getwd() #"C:/Users/Hyunsoo Kim/Documents/lecture/regression_analysis"[1] "C:/Users/Hyunsoo Kim/Desktop/senior_grade/blog/my-quarto-website"Regression Analysis
Basic
Code
R
0316 - 선형모형
예제를 통한 회귀분석
2장
표2.5
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")
dim(data_2.5) #14 2[1] 14 2head(data_2.5) Minutes Units
1 23 1
2 29 2
3 49 3
4 64 4
5 74 4
6 87 5X<-data_2.5$Units
Y<-data_2.5$Minutes표2.6
df<-data.frame(
#1:length(X),
Y,
X,
Y-mean(Y),
X-mean(X),
(Y-mean(Y))^2,
(X-mean(X))^2,
(Y-mean(Y))*(X-mean(X))
)
df Y X Y...mean.Y. X...mean.X. X.Y...mean.Y...2 X.X...mean.X...2
1 23 1 -74.2142857 -5 5.507760e+03 25
2 29 2 -68.2142857 -4 4.653189e+03 16
3 49 3 -48.2142857 -3 2.324617e+03 9
4 64 4 -33.2142857 -2 1.103189e+03 4
5 74 4 -23.2142857 -2 5.389031e+02 4
6 87 5 -10.2142857 -1 1.043316e+02 1
7 96 6 -1.2142857 0 1.474490e+00 0
8 97 6 -0.2142857 0 4.591837e-02 0
9 109 7 11.7857143 1 1.389031e+02 1
10 119 8 21.7857143 2 4.746173e+02 4
11 149 9 51.7857143 3 2.681760e+03 9
12 145 9 47.7857143 3 2.283474e+03 9
13 154 10 56.7857143 4 3.224617e+03 16
14 166 10 68.7857143 4 4.731474e+03 16
X.Y...mean.Y......X...mean.X..
1 371.07143
2 272.85714
3 144.64286
4 66.42857
5 46.42857
6 10.21429
7 0.00000
8 0.00000
9 11.78571
10 43.57143
11 155.35714
12 143.35714
13 227.14286
14 275.14286공분산(Covariance) - 식2.2
COV_XY<-sum((Y-mean(Y))*(X-mean(X))) / (length(X)-1) #136
### cov() 함수
cov(X,Y) #136[1] 136### 상관계수(correalationship)
### cor() 함수
cor(X,Y) #0.9936987 [1] 0.99369870321 - 선형모형
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")
x<-data_2.5$Units
y<-data_2.5$Minutes식 2.6
cor_xy<- COV_XY / (sd(x)*sd(y))
cor_xy[1] 0.9936987### cor() 함수
cor(x,y)[1] 0.9936987cor(y,x)[1] 0.9936987data_2.5 Minutes Units
1 23 1
2 29 2
3 49 3
4 64 4
5 74 4
6 87 5
7 96 6
8 97 6
9 109 7
10 119 8
11 149 9
12 145 9
13 154 10
14 166 10cor(data_2.5) Minutes Units
Minutes 1.0000000 0.9936987
Units 0.9936987 1.0000000그림 2.4
class(X)[1] "numeric"class(Y) #both numeric[1] "numeric"plot(X,Y, pch=19,xlab="Units",ylab="Minutes") 
표 2.3
data_2.3<-read.table("All_Data/p029a.txt",header=TRUE,sep="\t")
data_2.3 Y X
1 1 -7
2 14 -6
3 25 -5
4 34 -4
5 41 -3
6 46 -2
7 49 -1
8 50 0
9 49 1
10 46 2
11 41 3
12 34 4
13 25 5
14 14 6
15 1 7X<-data_2.3$X
Y<-data_2.3$Y그림 2.2
plot(X,Y)
cor(X,Y) # 0 완벽하게 2차함수의 형태도 0이 나옴(직선의 형태가 아닌것만)[1] 0표 2.4
data_2.4<-read.table("All_Data/p029b.txt",header=TRUE,sep="\t")그림 2.3
plot(data_2.4$X1,data_2.4$Y1, pch=19); abline(3,0.5) #기울기 3 절편0.5인 선을 추가해라 
plot(data_2.4$X2,data_2.4$Y2, pch=19); abline(3,0.5)
plot(data_2.4$X3,data_2.4$Y3, pch=19); abline(3,0.5)
plot(data_2.4$X4,data_2.4$Y4, pch=19); abline(3,0.5)
m<-matrix(1:4,ncol=2,byrow=TRUE) #2행의 매트릭스 생성
m [,1] [,2]
[1,] 1 2
[2,] 3 4layout(m)
plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)
#y~x -> y=ax+b 이러한 형태를 가지는 모형식이라는 의미
plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)
plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)
plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)
layout(1) #변환을 다시 하지 않으면 설정한 매트릭스의 비율로 그래프가 그려짐 해제 필요
# cor()
cor(data_2.4$X1,data_2.4$Y1) #0.8164205[1] 0.8164205cor(data_2.4$X2,data_2.4$Y2) #0.8162365[1] 0.8162365cor(data_2.4$X3,data_2.4$Y3) #0.8162867[1] 0.8162867cor(data_2.4$X4,data_2.4$Y4) #0.8165214[1] 0.8165214cor(data_2.4) #이렇게 한번에 할 수 있으나 가독성 떨어짐 Y1 X1 Y2 X2 Y3 X3 Y4
Y1 1.0000000 0.8164205 0.7500054 0.8164205 0.4687167 0.8164205 -0.4891162
X1 0.8164205 1.0000000 0.8162365 1.0000000 0.8162867 1.0000000 -0.3140467
Y2 0.7500054 0.8162365 1.0000000 0.8162365 0.5879193 0.8162365 -0.4780949
X2 0.8164205 1.0000000 0.8162365 1.0000000 0.8162867 1.0000000 -0.3140467
Y3 0.4687167 0.8162867 0.5879193 0.8162867 1.0000000 0.8162867 -0.1554718
X3 0.8164205 1.0000000 0.8162365 1.0000000 0.8162867 1.0000000 -0.3140467
Y4 -0.4891162 -0.3140467 -0.4780949 -0.3140467 -0.1554718 -0.3140467 1.0000000
X4 -0.5290927 -0.5000000 -0.7184365 -0.5000000 -0.3446610 -0.5000000 0.8165214
X4
Y1 -0.5290927
X1 -0.5000000
Y2 -0.7184365
X2 -0.5000000
Y3 -0.3446610
X3 -0.5000000
Y4 0.8165214
X4 1.0000000단순선형회귀모형 2.4
2.5 모수에 대한 추정
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")
x<-data_2.5$Units
y<-data_2.5$Minutes식 2.14 & 2.15
sum((y-mean(y))*(x-mean(x))) #1768[1] 1768sum((x-mean(x))^2) #114[1] 114beta1_hat<-sum((y-mean(y))*(x-mean(x))) / sum((x-mean(x))^2)
beta1_hat #15.50877[1] 15.50877beta0_hat <- mean(y) - (beta1_hat*mean(x))
beta0_hat #4.161654[1] 4.161654### 최소제곱회귀 방정식
# Minutes = 4.161654 + 15.50877 * Units
plot(beta0_hat+beta1_hat*x,pch=19);
# 4개의 고장 난 부품을 수리하는데 걸리는 에측시간
4.161654 + 15.50877 * 4 #66.19673[1] 66.19673units<-4
beta0_hat + beta1_hat*units[1] 66.19674### 적합값(Fitted value)
y_hat<-beta0_hat + beta1_hat*(x)
### 최소 제곱 잔차(residual)
e<-y-y_hat
e #합이 0이라는 특징이 존재 [1] 3.3295739 -6.1791980 -1.6879699 -2.1967419 7.8032581 5.2944862
[7] -1.2142857 -0.2142857 -3.7230576 -9.2318296 5.2593985 1.2593985
[13] -5.2493734 6.7506266sum(e) #1.278977e-13 0에 근사한 추지가 나옴[1] 1.278977e-13표2.7
df_2.7<-data.frame(
x=x,
y=y,
y_hat,
e
)
df_2.7 x y y_hat e
1 1 23 19.67043 3.3295739
2 2 29 35.17920 -6.1791980
3 3 49 50.68797 -1.6879699
4 4 64 66.19674 -2.1967419
5 4 74 66.19674 7.8032581
6 5 87 81.70551 5.2944862
7 6 96 97.21429 -1.2142857
8 6 97 97.21429 -0.2142857
9 7 109 112.72306 -3.7230576
10 8 119 128.23183 -9.2318296
11 9 149 143.74060 5.2593985
12 9 145 143.74060 1.2593985
13 10 154 159.24937 -5.2493734
14 10 166 159.24937 6.7506266### lm() 함수 (linear model)
# Minutes = beta0 + beta1 * Units + epsilon
# 모형식 : y~x
lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
4.162 15.509 res_lm<-lm(Minutes~Units,data=data_2.5)
res_lm
Call:
lm(formula = Minutes ~ Units, data = data_2.5)
Coefficients:
(Intercept) Units
4.162 15.509 # 리스트의 이름
names(res_lm) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model" # 회귀계수
res_lm$coefficients(Intercept) Units
4.161654 15.508772 coef(res_lm)(Intercept) Units
4.161654 15.508772 # 적합값
res_lm$fitted.values 1 2 3 4 5 6 7 8
19.67043 35.17920 50.68797 66.19674 66.19674 81.70551 97.21429 97.21429
9 10 11 12 13 14
112.72306 128.23183 143.74060 143.74060 159.24937 159.24937 fitted(res_lm) 1 2 3 4 5 6 7 8
19.67043 35.17920 50.68797 66.19674 66.19674 81.70551 97.21429 97.21429
9 10 11 12 13 14
112.72306 128.23183 143.74060 143.74060 159.24937 159.24937 # 최소제곱잔차
res_lm$residuals 1 2 3 4 5 6 7
3.3295739 -6.1791980 -1.6879699 -2.1967419 7.8032581 5.2944862 -1.2142857
8 9 10 11 12 13 14
-0.2142857 -3.7230576 -9.2318296 5.2593985 1.2593985 -5.2493734 6.7506266 resid(res_lm) 1 2 3 4 5 6 7
3.3295739 -6.1791980 -1.6879699 -2.1967419 7.8032581 5.2944862 -1.2142857
8 9 10 11 12 13 14
-0.2142857 -3.7230576 -9.2318296 5.2593985 1.2593985 -5.2493734 6.7506266 residuals(res_lm) 1 2 3 4 5 6 7
3.3295739 -6.1791980 -1.6879699 -2.1967419 7.8032581 5.2944862 -1.2142857
8 9 10 11 12 13 14
-0.2142857 -3.7230576 -9.2318296 5.2593985 1.2593985 -5.2493734 6.7506266 그림 2.5 - 산점도 + 회귀선
plot(beta0_hat+beta1_hat*x,pch=19);
#abline(beta0_hat,beta1_hat)
abline(res_lm)
0323
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")
res_lm <- lm(Minutes ~ Units, data = data_2.5)
res_lm
Call:
lm(formula = Minutes ~ Units, data = data_2.5)
Coefficients:
(Intercept) Units
4.162 15.509 res_lm_summ<-summary(res_lm)
res_lm_summ #Pr(>|t|) - p-value
Call:
lm(formula = Minutes ~ Units, data = data_2.5)
Residuals:
Min 1Q Median 3Q Max
-9.2318 -3.3415 -0.7143 4.7769 7.8033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.162 3.355 1.24 0.239
Units 15.509 0.505 30.71 8.92e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.392 on 12 degrees of freedom
Multiple R-squared: 0.9874, Adjusted R-squared: 0.9864
F-statistic: 943.2 on 1 and 12 DF, p-value: 8.916e-13# unit은 시간에 영향을준다 약15.5분 만큼씩
# coefficient에서 p-value에 대해서 알 수 있음
# beta_0는 0이라고 보면되느냐? p-value가 0.05보다 크기에 0328
2.7 신뢰구간
confint(res_lm) # beta_0,1의 95% 신뢰구간을 뽑아줌 2.5 % 97.5 %
(Intercept) -3.148482 11.47179
Units 14.408512 16.60903?confint #level = 1-alphastarting httpd help server ... doneconfint(res_lm, level=0.90) # 90%의 신뢰구간 5 % 95 %
(Intercept) -1.81810 10.14141
Units 14.60875 16.408792.8 예측
# 4개의 고장난 부품을 수리하는 데에 걸리는 시간 예측
x<-4
4.161654 + 15.508772 *4[1] 66.19674res_lm$coefficients[1]+(res_lm$coefficients[2]*x)(Intercept)
66.19674 ### predict()
df<-data.frame(Units=4)
predict(res_lm,newdata=df) # res_lm을 만들때 사용한 데이터형식으로 만들어주어야함 1
66.19674 res_lm_pred<-predict(res_lm,newdata=df,se.fit=TRUE)
### 예측값
res_lm_pred$fit 1
66.19674 ### 표준오차
res_lm_pred$se.fit # 평균반응에 대한 표준오차 [1] 1.759688### 예측한계
df<-data.frame(Units=4) #예제서는 4대기준
df<-data.frame(Units=1:10) #다른것도 보고 싶은경우
res_lm_pred_int_p<-predict(res_lm,newdata=df,interval="prediction")
### 신뢰한계
df<-data.frame(Units=1:10) #다른것도 보고 싶은경우
res_lm_pred_int_c<-predict(res_lm,newdata=df,interval="confidence") #둘의 차이를 보면 예측한계의 범위가 더큼
### 예측한계 & 신뢰한계
# 신뢰한계는 평균에서 멀어지만 오차의범위가 커지고 평균에 다가갈수록 오차가 줄어듬
plot(Minutes~Units,data=data_2.5,pch=19)
abline(res_lm,col="red",lwd=2)
lines(1:10,res_lm_pred_int_p[,"lwr"],col="darkgreen")
lines(1:10,res_lm_pred_int_p[,"upr"],col="darkgreen")
lines(1:10,res_lm_pred_int_c[,"lwr"],col="blue")
lines(1:10,res_lm_pred_int_c[,"upr"],col="blue")
2.9 적합성의 측정
res_lm_summ<-summary(res_lm)
res_lm_summ #Multiple R-squared:0.9874 -> 반응변수의 전체변이중 98.94%가 예측변수에 의해 설명된다
Call:
lm(formula = Minutes ~ Units, data = data_2.5)
Residuals:
Min 1Q Median 3Q Max
-9.2318 -3.3415 -0.7143 4.7769 7.8033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.162 3.355 1.24 0.239
Units 15.509 0.505 30.71 8.92e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.392 on 12 degrees of freedom
Multiple R-squared: 0.9874, Adjusted R-squared: 0.9864
F-statistic: 943.2 on 1 and 12 DF, p-value: 8.916e-13# 만약 R-squared가 1이면 완벽한 선형의 관계 100%라는 것을 의미한다.
# R-squared는 변수가 들어갈수록 커지기에 adjust R-squared를 사용 추후 설명 2.10 원점을 통과하는 회귀선
# Minutes = beta1 + Units + epsilon
res_lm_no<-lm(Minutes~Units-1,data=data_2.5)
res_lm_no
Call:
lm(formula = Minutes ~ Units - 1, data = data_2.5)
Coefficients:
Units
16.07 summary(res_lm_no)
Call:
lm(formula = Minutes ~ Units - 1, data = data_2.5)
Residuals:
Min 1Q Median 3Q Max
-9.5955 -2.4733 0.4417 5.0243 9.7023
Coefficients:
Estimate Std. Error t value Pr(>|t|)
Units 16.0744 0.2213 72.62 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.502 on 13 degrees of freedom
Multiple R-squared: 0.9975, Adjusted R-squared: 0.9974
F-statistic: 5274 on 1 and 13 DF, p-value: < 2.2e-16coef(summary(res_lm_no)) #rsquared=0.9975 Estimate Std. Error t value Pr(>|t|)
Units 16.07443 0.2213341 72.62519 2.380325e-182.11
y<-rnorm(30)
t.test(y,mu=0)
One Sample t-test
data: y
t = -0.67175, df = 29, p-value = 0.5071
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-0.5134310 0.2595474
sample estimates:
mean of x
-0.1269418 summary(lm(y~1))
Call:
lm(formula = y ~ 1)
Residuals:
Min 1Q Median 3Q Max
-1.96463 -0.67286 0.02797 0.69205 2.12145
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1269 0.1890 -0.672 0.507
Residual standard error: 1.035 on 29 degrees of freedom3장 다중선형회귀
3.3 사례 감독자 직무수행능력 데이터
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
dim(data_3.3)[1] 30 7class(data_3.3)[1] "data.frame"sapply(data_3.3,class) #all numeric Y X1 X2 X3 X4 X5 X6
"numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" summary(data_3.3) #모든변수가 numeric이면 분위수도 보여준다 Y X1 X2 X3 X4
Min. :40.00 Min. :37.0 Min. :30.00 Min. :34.00 Min. :43.00
1st Qu.:58.75 1st Qu.:58.5 1st Qu.:45.00 1st Qu.:47.00 1st Qu.:58.25
Median :65.50 Median :65.0 Median :51.50 Median :56.50 Median :63.50
Mean :64.63 Mean :66.6 Mean :53.13 Mean :56.37 Mean :64.63
3rd Qu.:71.75 3rd Qu.:77.0 3rd Qu.:62.50 3rd Qu.:66.75 3rd Qu.:71.00
Max. :85.00 Max. :90.0 Max. :83.00 Max. :75.00 Max. :88.00
X5 X6
Min. :49.00 Min. :25.00
1st Qu.:69.25 1st Qu.:35.00
Median :77.50 Median :41.00
Mean :74.77 Mean :42.93
3rd Qu.:80.00 3rd Qu.:47.75
Max. :92.00 Max. :72.00 ### 산점도 행렬
plot(data_3.3)
3.4 모수 추정
lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)
Call:
lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data_3.3)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
10.78708 0.61319 -0.07305 0.32033 0.08173 0.03838
X6
-0.21706 lm(Y~.,data=data_3.3) # X1+X2+X3+X4+X5+X6쓰는 것이 아니라 .을 써서 모든 변수를 다써줌
Call:
lm(formula = Y ~ ., data = data_3.3)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
10.78708 0.61319 -0.07305 0.32033 0.08173 0.03838
X6
-0.21706 3.5 회귀계수에 대한 해석
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
lm(Y~X1+X2,data=data_3.3)
Call:
lm(formula = Y ~ X1 + X2, data = data_3.3)
Coefficients:
(Intercept) X1 X2
15.32762 0.78034 -0.05016 # (Intercept) X1 X2
# 15.32762 0.78034 -0.05016
# 1) Y에서 X1 효과 제거
m1<-lm(Y~X1,data=data_3.3) # y prime
m1$residuals # x1이 설명하지 못한값 / x1의 효과를 제거한 값 1 2 3 4 5 6
-9.86142016 0.32865220 3.80099328 -0.91673799 7.76411473 -12.87985944
7 8 9 10 11 12
-6.93517726 0.02794419 -4.25432454 6.59248165 9.62936020 7.34709147
13 14 15 16 17 18
7.83787183 -9.00893436 4.51872455 -1.29120309 -4.51815400 5.34709147
19 20 21 22 23 24
-2.19900672 -8.14368889 5.43928784 3.59248165 -11.18056744 -2.29688270
25 26 27 28 29 30
7.87475038 -6.48127545 7.02794419 -9.38907907 6.48184600 5.74567546 # 2) X2에서 X1 효과 제거
m2<-lm(X2~X1,data=data_3.3) # x2 prime
m2$residuals 1 2 3 4 5 6
-15.1300345 -0.7994502 13.1223579 -6.2864182 -2.9818979 1.8178376
7 8 9 10 11 12
-11.3385461 -7.4428019 10.9659742 -5.2603543 6.8439016 -2.7473223
13 14 15 16 17 18
6.2266138 21.4529422 -4.4688659 -15.1382816 1.4268783 15.2526777
19 20 21 22 23 24
-8.8776421 19.2787417 -6.4866827 1.7396457 -0.8255141 4.0524132
25 26 27 28 29 30
-4.6691304 7.5311341 0.5571981 -4.2082264 8.4268783 -22.0340258 # 3) X1의 효과가 제거된 Y와 X2의 적합 - 원점을 지나는 회귀선
lm(m1$residuals~m2$residuals-1) # 원점을 지나면 -1를 하고 진행 // -3.25e-17
Call:
lm(formula = m1$residuals ~ m2$residuals - 1)
Coefficients:
m2$residuals
-0.05016 # 다른 효과 없이(다른값이 고정) Y에 영향을 주는 순수한 X2의 값
# m2$residuals : -0.05016 == X2 : -0.05016
### 단위길이 척도화 - 잘사용하지않음
fn_scaling_len<-function(x){
x0<-x-mean(x)
x0/sqrt(sum(x0^2))
}
data_3.3_len<-sapply(data_3.3, fn_scaling_len)
data_3.3_len<-data.frame(data_3.3_len)
summary(data_3.3_len) Y X1 X2 X3
Min. :-0.37579 Min. :-0.41282 Min. :-0.35109 Min. :-0.353871
1st Qu.:-0.08975 1st Qu.:-0.11297 1st Qu.:-0.12344 1st Qu.:-0.148193
Median : 0.01322 Median :-0.02231 Median :-0.02479 Median : 0.002109
Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.000000
3rd Qu.: 0.10857 3rd Qu.: 0.14504 3rd Qu.: 0.14216 3rd Qu.: 0.164278
Max. : 0.31070 Max. : 0.32635 Max. : 0.45328 Max. : 0.294804
X4 X5 X6
Min. :-0.38637 Min. :-0.48356 Min. :-0.32367
1st Qu.:-0.11401 1st Qu.:-0.10353 1st Qu.:-0.14318
Median :-0.02024 Median : 0.05130 Median :-0.03489
Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
3rd Qu.: 0.11371 3rd Qu.: 0.09821 3rd Qu.: 0.08693
Max. : 0.41733 Max. : 0.32341 Max. : 0.52461 lm(Y~.,data=data_3.3_len)
Call:
lm(formula = Y ~ ., data = data_3.3_len)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
-1.259e-16 6.707e-01 -7.343e-02 3.089e-01 6.981e-02 3.120e-02
X6
-1.835e-01 ### 표준화
# scale()
data_3.3_std<-scale(data_3.3)
#summary(data_3.3_std)
#sapply(data_3.3_std, sd, na.rm=T)
#class(data_3.3_std) #"matrix"
data_3.3_std<-data.frame(data_3.3_std)
#class(data_3.3_std) #"data.frame"
#sapply(data_3.3_std, sd, na.rm=T)
lm(Y~.,data=data_3.3_std) # beta게수 구하기
Call:
lm(formula = Y ~ ., data = data_3.3_std)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
-7.717e-16 6.707e-01 -7.343e-02 3.089e-01 6.981e-02 3.120e-02
X6
-1.835e-01 3.7 최소제곱추정량의 성질
res_lm<-lm(Y~.,data=data_3.3)
res_lm
Call:
lm(formula = Y ~ ., data = data_3.3)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
10.78708 0.61319 -0.07305 0.32033 0.08173 0.03838
X6
-0.21706 summary(res_lm)
Call:
lm(formula = Y ~ ., data = data_3.3)
Residuals:
Min 1Q Median 3Q Max
-10.9418 -4.3555 0.3158 5.5425 11.5990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.78708 11.58926 0.931 0.361634
X1 0.61319 0.16098 3.809 0.000903 ***
X2 -0.07305 0.13572 -0.538 0.595594
X3 0.32033 0.16852 1.901 0.069925 .
X4 0.08173 0.22148 0.369 0.715480
X5 0.03838 0.14700 0.261 0.796334
X6 -0.21706 0.17821 -1.218 0.235577
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05m1<-summary(lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)) #Adjusted R-squared: 0.6628
m2<-summary(lm(Y~X1+X2+X3+X4+X5,data=data_3.3)) #Adjusted R-squared: 0.6561
# X6가 들어가는 것이 더 좋은 모델
m1$adj.r.squared[1] 0.662846m2$adj.r.squared #summary에서 보다 더 정확하게 수치가 나옴 [1] 0.65605393.9 개별 회귀계수들에 대한 추론
res_lm_summ<-summary(res_lm)
res_lm_summ #p-value의 존재는 무언가를 검정했다라는 반증
Call:
lm(formula = Y ~ ., data = data_3.3)
Residuals:
Min 1Q Median 3Q Max
-10.9418 -4.3555 0.3158 5.5425 11.5990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.78708 11.58926 0.931 0.361634
X1 0.61319 0.16098 3.809 0.000903 ***
X2 -0.07305 0.13572 -0.538 0.595594
X3 0.32033 0.16852 1.901 0.069925 .
X4 0.08173 0.22148 0.369 0.715480
X5 0.03838 0.14700 0.261 0.796334
X6 -0.21706 0.17821 -1.218 0.235577
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05# p-value<0.05 H_1 귀무가설 채택
# p-value>0.05 H_0 영가설 채택 // X1을 제외하고는 영가설 유의한 의미가 없음(Y에영향주는)
# 모두다 0이라는 가설을 가지고 분모 분자의 오차가 카이제곱을 따르고 거기서 나온 통계량
# F-분포 자유도는 분자 분모 두개를 가짐 //모아서 계산을 하기에 각각 계산하는것과 결과다름
# 영가설-모든 회귀계수가 0이다.
# 대립가설-적어도 하나는 0이 아니다. p-value: 1.24e-05 <0.05 대립가설 채택
# p-value가 0.05보다 작으면 대립가설 채택!!!!!! 기억해
# 회귀계수에 대한 신뢰구간 - 95% 신뢰한계
confint(res_lm) #-13.18712881 ~ 34.7612816 2.5 % 97.5 %
(Intercept) -13.18712881 34.7612816
X1 0.28016866 0.9462066
X2 -0.35381806 0.2077178
X3 -0.02827872 0.6689430
X4 -0.37642935 0.5398936
X5 -0.26570179 0.3424647
X6 -0.58571106 0.1515977#X1 0.28016866 0.9462066 사이에 0이 들어가있으면 영향을 준다라느걸 의미
#X2 -0.35381806 0.2077178 p-value없이도 알 수 있음
#X5가 가장 영향이 적음 p-value가 가장 크기에(영향 효과의 크기를 비교할때)
#p-value가 작을 수록 영향을 많이 준다 beta값을 보는 것이 아닌 p-value를 보는 것 중요
#가장 의미있는 변수?->p-value가장 작은거 // 대립가설채택 Y에 영향을 가장3.10 선형모형에서의 가설검정
3.10.1 모든 회귀예수들이 0인가에 대한 검정
# H_0: beta_1:beta_6=0
model_reduce<-lm(Y~1,data=data_3.3)
model_full<-lm(Y~.,data=data_3.3)
anova(model_reduce,model_full)Analysis of Variance Table
Model 1: Y ~ 1
Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6
Res.Df RSS Df Sum of Sq F Pr(>F)
1 29 4297
2 23 1149 6 3148 10.502 1.24e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#대립가설 = 완전모형이 적절하다 / 1.24e-05 *** < 0.05
#의미 있는 예측 변수가 한개 이상 존재한다
summary(model_full) #summary에서 beta_1~beta_6까지 모두가 0이라는 가설로 진행을 이미함
Call:
lm(formula = Y ~ ., data = data_3.3)
Residuals:
Min 1Q Median 3Q Max
-10.9418 -4.3555 0.3158 5.5425 11.5990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.78708 11.58926 0.931 0.361634
X1 0.61319 0.16098 3.809 0.000903 ***
X2 -0.07305 0.13572 -0.538 0.595594
X3 0.32033 0.16852 1.901 0.069925 .
X4 0.08173 0.22148 0.369 0.715480
X5 0.03838 0.14700 0.261 0.796334
X6 -0.21706 0.17821 -1.218 0.235577
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05#F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05
#예상// 가장의미있는변수? -> X1 이유?-> p-value 0.000903 로 가장 작기에 영향많이 줄것으로 예측 3.10.2 회귀계수들의 부분집합이 0인가에 대한 검정
model_reduce<-lm(Y~X1+X3,data=data_3.3)
model_full<-lm(Y~.,data=data_3.3)
anova(model_reduce,model_full) #0.7158 > 0.05Analysis of Variance Table
Model 1: Y ~ X1 + X3
Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27 1254.7
2 23 1149.0 4 105.65 0.5287 0.7158#영가설은 H_0: b_2=b_4=b_5=b_6=0 이라는 사실을 알 수 있다
#b_1&b_3는 반응변수에 유의한 반응을 준다라는 것도 연계하여 알 수 있다 3.10.3 회귀계수들의 통일성에 대한 검정
#해당 조건이 주어지고 만족할 때 beta_1=beta_3은 맞는가?
model_reduce<-lm(Y~I(X1-X3),data=data_3.3) #I를 씌우면 새로운 변수를 만든것과 동일
# X1-X3를 한 그자체를 분석하라는 의미//본래는 X1-X3 해서 새로운 변수를 만들어서 해야함
model_full<-lm(Y~X1+X3,data=data_3.3)
anova(model_reduce,model_full) Analysis of Variance Table
Model 1: Y ~ I(X1 - X3)
Model 2: Y ~ X1 + X3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 3846.7
2 27 1254.6 1 2592 55.78 4.925e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#install.packages("car")
library(car)Warning: package 'car' was built under R version 4.3.2Loading required package: carDataWarning: package 'carData' was built under R version 4.3.2model_full<-lm(Y~X1+X3,data=data_3.3)
car::linearHypothesis(model_full,c("X1=X3"))Linear hypothesis test
Hypothesis:
X1 - X3 = 0
Model 1: restricted model
Model 2: Y ~ X1 + X3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 1424.6
2 27 1254.7 1 169.95 3.6572 0.06649 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.10.4 제약조건하에서 회귀계수에 대한 추정과 검정
# H_0: beta_1+beta_3=1 | beta_2=beta_3:beta_6=0
model_full<-lm(Y~X1+X3,data=data_3.3)
car::linearHypothesis(model_full,c("X1 + X3 = 1"))Linear hypothesis test
Hypothesis:
X1 + X3 = 1
Model 1: restricted model
Model 2: Y ~ X1 + X3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 1329.5
2 27 1254.7 1 74.898 1.6118 0.2151# x1의 효과가 증가하면 x3의 효과는 감소한다 상대적인 관계 (반대로도 가능)3.11 예측
model_full<-lm(Y~.,data_3.3)
# 예측값 - 적합값
model_full$fitted.values 1 2 3 4 5 6 7 8
51.11030 61.35277 69.93944 61.22684 74.45380 53.94185 67.14841 70.09701
9 10 11 12 13 14 15 16
79.53099 59.19846 57.92572 55.40103 59.58168 70.21401 76.54933 84.54785
17 18 19 20 21 22 23 24
76.15013 61.39736 68.01656 55.62014 42.60324 63.81902 63.66400 44.62475
25 26 27 28 29 30
57.31710 67.84347 75.14036 56.04535 77.66053 76.87850 # 예측한계(Prediction Limits)
predict(model_full,newdata = data_3.3,interval = "prediction") fit lwr upr
1 51.11030 34.16999 68.05060
2 61.35277 46.34536 76.36018
3 69.93944 53.94267 85.93622
4 61.22684 45.44586 77.00783
5 74.45380 59.17630 89.73129
6 53.94185 37.21813 70.66557
7 67.14841 51.64493 82.65189
8 70.09701 54.54384 85.65017
9 79.53099 62.71383 96.34814
10 59.19846 44.03506 74.36185
11 57.92572 42.00674 73.84470
12 55.40103 39.79333 71.00873
13 59.58168 43.39853 75.76483
14 70.21401 52.06636 88.36167
15 76.54933 60.79444 92.30422
16 84.54785 66.41374 102.68197
17 76.15013 59.99991 92.30036
18 61.39736 43.23384 79.56088
19 68.01656 52.44673 83.58639
20 55.62014 39.63744 71.60284
21 42.60324 26.35046 58.85603
22 63.81902 48.40145 79.23659
23 63.66400 48.56222 78.76578
24 44.62475 27.25435 61.99514
25 57.31710 41.29380 73.34041
26 67.84347 49.98605 85.70089
27 75.14036 59.31975 90.96097
28 56.04535 40.18723 71.90348
29 77.66053 61.97564 93.34541
30 76.87850 60.27441 93.48258# 신뢰한계(Confidence limits)
predict(model_full,newdata = data_3.3,interval = "confidence") fit lwr upr
1 51.11030 42.55502 59.66557
2 61.35277 57.97029 64.73524
3 69.93944 63.44979 76.42909
4 61.22684 55.28897 67.16471
5 74.45380 70.02428 78.88332
6 53.94185 45.82386 62.05984
7 67.14841 61.99316 72.30367
8 70.09701 64.79421 75.39980
9 79.53099 71.22222 87.83975
10 59.19846 55.18008 63.21683
11 57.92572 51.63028 64.22116
12 55.40103 49.94035 60.86171
13 59.58168 52.64531 66.51805
14 70.21401 59.46431 80.96372
15 76.54933 70.68118 82.41748
16 84.54785 73.82102 95.27468
17 76.15013 69.29093 83.00933
18 61.39736 50.62090 72.17383
19 68.01656 62.66507 73.36805
20 55.62014 49.16527 62.07502
21 42.60324 35.50593 49.70055
22 63.81902 58.92819 68.70985
23 63.66400 59.88479 67.44321
24 44.62475 35.24662 54.00288
25 57.31710 50.76233 63.87187
26 67.84347 57.59134 78.09561
27 75.14036 69.09798 81.18275
28 56.04535 49.90540 62.18531
29 77.66053 71.98300 83.33806
30 76.87850 69.00992 84.74707부록 : 행렬을 이용한 회귀계수 추정
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
Y<-data_3.3$Y
X<-data_3.3[,-1]
X<-cbind(1,X)
X<-as.matrix(X)
#beta_hat<-solve(t(X) %*% X) %*% t(X) %*% Y # %*%행렬 계산
P<-solve(t(X) %*% X) %*% t(X)
beta_hat<- P %*% Y
lm(Y~.,data_3.3)
Call:
lm(formula = Y ~ ., data = data_3.3)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
10.78708 0.61319 -0.07305 0.32033 0.08173 0.03838
X6
-0.21706 4장 회귀진단: 모형위반의 검출
4.3 다양한 유형의 잔차들
# 표준화 잔차
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
res_lm<-lm(Y~.,data=data_3.3)
class(res_lm)[1] "lm"mode(res_lm)[1] "list"names(res_lm) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model" res_lm$fitted.values 1 2 3 4 5 6 7 8
51.11030 61.35277 69.93944 61.22684 74.45380 53.94185 67.14841 70.09701
9 10 11 12 13 14 15 16
79.53099 59.19846 57.92572 55.40103 59.58168 70.21401 76.54933 84.54785
17 18 19 20 21 22 23 24
76.15013 61.39736 68.01656 55.62014 42.60324 63.81902 63.66400 44.62475
25 26 27 28 29 30
57.31710 67.84347 75.14036 56.04535 77.66053 76.87850 str(res_lm)List of 12
$ coefficients : Named num [1:7] 10.7871 0.6132 -0.0731 0.3203 0.0817 ...
..- attr(*, "names")= chr [1:7] "(Intercept)" "X1" "X2" "X3" ...
$ residuals : Named num [1:30] -8.11 1.647 1.061 -0.227 6.546 ...
..- attr(*, "names")= chr [1:30] "1" "2" "3" "4" ...
$ effects : Named num [1:30] -354.011 54.107 2.742 11.715 -0.971 ...
..- attr(*, "names")= chr [1:30] "(Intercept)" "X1" "X2" "X3" ...
$ rank : int 7
$ fitted.values: Named num [1:30] 51.1 61.4 69.9 61.2 74.5 ...
..- attr(*, "names")= chr [1:30] "1" "2" "3" "4" ...
$ assign : int [1:7] 0 1 2 3 4 5 6
$ qr :List of 5
..$ qr : num [1:30, 1:7] -5.477 0.183 0.183 0.183 0.183 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:7] "(Intercept)" "X1" "X2" "X3" ...
.. ..- attr(*, "assign")= int [1:7] 0 1 2 3 4 5 6
..$ qraux: num [1:7] 1.18 1 1.29 1.1 1.07 ...
..$ pivot: int [1:7] 1 2 3 4 5 6 7
..$ tol : num 1e-07
..$ rank : int 7
..- attr(*, "class")= chr "qr"
$ df.residual : int 23
$ xlevels : Named list()
$ call : language lm(formula = Y ~ ., data = data_3.3)
$ terms :Classes 'terms', 'formula' language Y ~ X1 + X2 + X3 + X4 + X5 + X6
.. ..- attr(*, "variables")= language list(Y, X1, X2, X3, X4, X5, X6)
.. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:7] "Y" "X1" "X2" "X3" ...
.. .. .. ..$ : chr [1:6] "X1" "X2" "X3" "X4" ...
.. ..- attr(*, "term.labels")= chr [1:6] "X1" "X2" "X3" "X4" ...
.. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(Y, X1, X2, X3, X4, X5, X6)
.. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:7] "Y" "X1" "X2" "X3" ...
$ model :'data.frame': 30 obs. of 7 variables:
..$ Y : num [1:30] 43 63 71 61 81 43 58 71 72 67 ...
..$ X1: num [1:30] 51 64 70 63 78 55 67 75 82 61 ...
..$ X2: num [1:30] 30 51 68 45 56 49 42 50 72 45 ...
..$ X3: num [1:30] 39 54 69 47 66 44 56 55 67 47 ...
..$ X4: num [1:30] 61 63 76 54 71 54 66 70 71 62 ...
..$ X5: num [1:30] 92 73 86 84 83 49 68 66 83 80 ...
..$ X6: num [1:30] 45 47 48 35 47 34 35 41 31 41 ...
..- attr(*, "terms")=Classes 'terms', 'formula' language Y ~ X1 + X2 + X3 + X4 + X5 + X6
.. .. ..- attr(*, "variables")= language list(Y, X1, X2, X3, X4, X5, X6)
.. .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:7] "Y" "X1" "X2" "X3" ...
.. .. .. .. ..$ : chr [1:6] "X1" "X2" "X3" "X4" ...
.. .. ..- attr(*, "term.labels")= chr [1:6] "X1" "X2" "X3" "X4" ...
.. .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. .. ..- attr(*, "predvars")= language list(Y, X1, X2, X3, X4, X5, X6)
.. .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
.. .. .. ..- attr(*, "names")= chr [1:7] "Y" "X1" "X2" "X3" ...
- attr(*, "class")= chr "lm"# 잔차
res_lm$residuals 1 2 3 4 5 6
-8.1102953 1.6472337 1.0605589 -0.2268416 6.5462010 -10.9418499
7 8 9 10 11 12
-9.1484140 0.9029929 -7.5309862 7.8015424 6.0742817 11.5989723
13 14 15 16 17 18
9.4183197 -2.2140147 0.4506705 -3.5478519 -2.1501319 3.6026355
19 20 21 22 23 24
-3.0165587 -5.6201442 7.3967582 0.1809831 -10.6639999 -4.6247464
25 26 27 28 29 30
5.6828983 -1.8434727 2.8596385 -8.0453540 7.3394730 5.1215016 resid(res_lm) #실제값에서 예측된 값을 뺸값 1 2 3 4 5 6
-8.1102953 1.6472337 1.0605589 -0.2268416 6.5462010 -10.9418499
7 8 9 10 11 12
-9.1484140 0.9029929 -7.5309862 7.8015424 6.0742817 11.5989723
13 14 15 16 17 18
9.4183197 -2.2140147 0.4506705 -3.5478519 -2.1501319 3.6026355
19 20 21 22 23 24
-3.0165587 -5.6201442 7.3967582 0.1809831 -10.6639999 -4.6247464
25 26 27 28 29 30
5.6828983 -1.8434727 2.8596385 -8.0453540 7.3394730 5.1215016 ### 내적 표준화잔차
rstandard(res_lm) 1 2 3 4 5 6
-1.41498026 0.23955370 0.16744867 -0.03512080 0.97184596 -1.86133876
7 8 9 10 11 12
-1.38317210 0.13709194 -1.29490454 1.14799070 0.95218982 1.76906521
13 14 15 16 17 18
1.51371017 -0.46212316 0.06961486 -0.73868563 -0.34446368 0.75418016
19 20 21 22 23 24
-0.45861365 -0.88618779 1.19699287 0.02717120 -1.56184734 -0.85286680
25 26 27 28 29 30
0.89948517 -0.36581416 0.44430497 -1.25422677 1.12683185 0.85971512 ### 외적 표준화잔차
MASS::studres(res_lm) 1 2 3 4 5 6
-1.44835328 0.23458097 0.16386794 -0.03434974 0.97062209 -1.97526518
7 8 9 10 11 12
-1.41280382 0.13413337 -1.31529351 1.15637546 0.95017640 1.86145176
13 14 15 16 17 18
1.56019127 -0.45407837 0.06809185 -0.73117411 -0.33776450 0.74689589
19 20 21 22 23 24
-0.45059801 -0.88189556 1.20894332 0.02657438 -1.61559196 -0.84763116
25 26 27 28 29 30
0.89560731 -0.35881868 0.43641573 -1.27088888 1.13380428 0.85466249 redsid_df<-data.frame(
Y=data_3.3$Y,
Y_hat=res_lm$fitted.values,
resid=resid(res_lm),
rstandard=rstandard(res_lm),
studres=MASS::studres(res_lm)
)
redsid_df Y Y_hat resid rstandard studres
1 43 51.11030 -8.1102953 -1.41498026 -1.44835328
2 63 61.35277 1.6472337 0.23955370 0.23458097
3 71 69.93944 1.0605589 0.16744867 0.16386794
4 61 61.22684 -0.2268416 -0.03512080 -0.03434974
5 81 74.45380 6.5462010 0.97184596 0.97062209
6 43 53.94185 -10.9418499 -1.86133876 -1.97526518
7 58 67.14841 -9.1484140 -1.38317210 -1.41280382
8 71 70.09701 0.9029929 0.13709194 0.13413337
9 72 79.53099 -7.5309862 -1.29490454 -1.31529351
10 67 59.19846 7.8015424 1.14799070 1.15637546
11 64 57.92572 6.0742817 0.95218982 0.95017640
12 67 55.40103 11.5989723 1.76906521 1.86145176
13 69 59.58168 9.4183197 1.51371017 1.56019127
14 68 70.21401 -2.2140147 -0.46212316 -0.45407837
15 77 76.54933 0.4506705 0.06961486 0.06809185
16 81 84.54785 -3.5478519 -0.73868563 -0.73117411
17 74 76.15013 -2.1501319 -0.34446368 -0.33776450
18 65 61.39736 3.6026355 0.75418016 0.74689589
19 65 68.01656 -3.0165587 -0.45861365 -0.45059801
20 50 55.62014 -5.6201442 -0.88618779 -0.88189556
21 50 42.60324 7.3967582 1.19699287 1.20894332
22 64 63.81902 0.1809831 0.02717120 0.02657438
23 53 63.66400 -10.6639999 -1.56184734 -1.61559196
24 40 44.62475 -4.6247464 -0.85286680 -0.84763116
25 63 57.31710 5.6828983 0.89948517 0.89560731
26 66 67.84347 -1.8434727 -0.36581416 -0.35881868
27 78 75.14036 2.8596385 0.44430497 0.43641573
28 48 56.04535 -8.0453540 -1.25422677 -1.27088888
29 85 77.66053 7.3394730 1.12683185 1.13380428
30 82 76.87850 5.1215016 0.85971512 0.854662494.5 모형을 적합깅 이전의 그래프
4.5.1 일차원 그래프
a<-rnorm(100,70,10) #연속형 데이터
# 히스토그램
hist(a)
hist(a,breaks=5) #범위를 조절 막대의 5번 자름 
# 줄기 잎 그림
stem(a)
The decimal point is 1 digit(s) to the right of the |
3 | 8
4 |
5 | 01146678999999
6 | 001122223444456666666777888899999
7 | 000111222333333445555667777778889
8 | 011122333334557789
9 | 2stem(round(a)) #줄기잎을 그릴때는 반올림을 하고 항상 진행
The decimal point is 1 digit(s) to the right of the |
3 | 8
4 |
5 | 01146678999999
6 | 001122223444456666666777888899999
7 | 000111222333333445555667777778889
8 | 011122333334557789
9 | 2stem(round(a),scale=2) #scale을 2배로 늘려라 5기준으로 반으로 잘라서
The decimal point is 1 digit(s) to the right of the |
3 | 8
4 |
4 |
5 | 0114
5 | 6678999999
6 | 0011222234444
6 | 56666666777888899999
7 | 00011122233333344
7 | 5555667777778889
8 | 011122333334
8 | 557789
9 | 2# 모든데이터를 볼 수있는 장점 데이터가 많으면 구림
# 점플롯
idx<-rep(1,length(a)) #a의 갯수에 맞춰서 1를 반복
plot(idx,a)
plot(jitter(idx),a,xlim=c(0.5,1.5))
# 상자그림
boxplot(a) #사분위수에 대해서 알 수 있음
# 상자그림 + 점플롯
boxplot(a)
points(jitter(idx),a)
4.5.2 이차원 그래프
data_4.1<-read.table("All_Data/p103.txt",header=T,sep="\t")
data_4.1 Y X1 X2
1 12.37 2.23 9.66
2 12.66 2.57 8.94
3 12.00 3.87 4.40
4 11.93 3.10 6.64
5 11.06 3.39 4.91
6 13.03 2.83 8.52
7 13.13 3.02 8.04
8 11.44 2.14 9.05
9 12.86 3.04 7.71
10 10.84 3.26 5.11
11 11.20 3.39 5.05
12 11.56 2.35 8.51
13 10.83 2.76 6.59
14 12.63 3.90 4.90
15 12.46 3.16 6.96class(data_4.1) #data.frame[1] "data.frame"# 산점도 행렬
plot(data_4.1)
cor(data_4.1) #상관계수 Y X1 X2
Y 1.000000000 0.002497966 0.4340688
X1 0.002497966 1.000000000 -0.8997765
X2 0.434068758 -0.899776481 1.0000000pairs(data_4.1)
# Correlation panel
panel.cor<-function(x,y){
par(usr=c(0,1,0,1))
r<-round(cor(x,y),digits = 3)
text(0.5,0.5,r,cex=1.5)
}
pairs(data_4.1,lower.panel = panel.cor)
# 회전도표, 동적 그래프(3차원)
#install.packages("rgl")
library(rgl)Warning: package 'rgl' was built under R version 4.3.2plot3d(x=data_4.1$X1,y=data_4.1$X2,z=data_4.1$Y) 4.7 선형성과 정규성 가정에 대한 검토
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
res_lm<-lm(Y~X1+X3,data=data_3.3)
layout(matrix(1:4,nrow=2,byrow=T))
plot(res_lm) #회귀진단 플랏이 나옴 
layout(1)
# 1. 표준화잔차의 정규확률분포
plot(res_lm,2)
# 2. 표준화잔차 대 각 예측변수들의 산점도
plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서 
# 3. 표준화잔차 대 적합값의 플롯
# 4. 표준화잔차의 인덱스 플롯0502 - 선형모형
library(dplyr)
Attaching package: 'dplyr'The following object is masked from 'package:car':
recodeThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union4.7 선형성과 정규성 가정에 대한 검토
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
res_lm<-lm(Y~X1+X3,data=data_3.3) #두개의 변수만 의미있다고 가정하고 진행
# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(res_lm) #회귀진단 플랏이 나옴 / 총6개임
layout(1) #다시 한개의 플랏만 그리도록
# 1. 표준화잔차의 정규확률분포
plot(res_lm,2) #QQ-plot y=x 기울기의 직선위에 점들이 있어야 한다 눈대중으로
# 2. 표준화잔차 대 각 예측변수들의 산점도
plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서
#랜덤하게 데이터가 흩어져 있어야 한다
# 내적 표준화잔차
plot(data_3.3$X1,rstandard(res_lm))
plot(data_3.3$X3,rstandard(res_lm)) #각각의 잔차들이 랜덤하게 잘 퍼져야함 
# 3. 표준화잔차 대 적합값의 플롯
plot(res_lm,1) #잔차와 적합값은 상관성이 없어야하며 랜덤하게 퍼져야함 
# 4. 표준화잔차의 인덱스 플롯
plot(res_lm,5)
data_2.4<-read.table("All_Data/p029b.txt",header=T,sep="\t")
m<-matrix(1:4,ncol=2,byrow=TRUE)
layout(m)
plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)
plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)
plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)
plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)
layout(1)Figure4.4
res_lm<-lm(Y1~X1,data=data_2.4)
#1번플랏은 적당히 잘퍼짐,2번플랏은 어느정도 선형성 있음(데이터적어서그런거임)
res_lm<-lm(Y2~X2,data=data_2.4)
res_lm<-lm(Y3~X3,data=data_2.4)
res_lm<-lm(Y4~X4,data=data_2.4)
res_lm
Call:
lm(formula = Y4 ~ X4, data = data_2.4)
Coefficients:
(Intercept) X4
3.0017 0.4999 layout(matrix(1:4,nrow=2,byrow=T))
plot(res_lm) Warning: not plotting observations with leverage one:
8
layout(1)4.8 지레점, 영향력, 특이값
# 사례: 뉴욕 강 데이터
# agr-농업, forest-숲, rsdntial-주거, comlndl-산업, nitrogen-질소
data_1.9<-read.table("All_Data/p010.txt",header=T,sep="\t")
head(data_1.9) River Agr Forest Rsdntial ComIndl Nitrogen
1 Olean 26 63 1.2 0.29 1.10
2 Cassadaga 29 57 0.7 0.09 1.01
3 Oatka 54 26 1.8 0.58 1.90
4 Neversink 2 84 1.9 1.98 1.00
5 Hackensack 3 27 29.4 3.11 1.99
6 Wappinger 19 61 3.4 0.56 1.42plot(data_1.9[-1],pch=19) #river라는 첫번째 컬럼을 제외하고 진행 
res_1<-lm(Nitrogen~.,data=data_1.9[-1]) #모든데이터 사용
res_2<-lm(Nitrogen~.,data=data_1.9[-4,-1]) #4번쨰 데이터 제외
res_3<-lm(Nitrogen~.,data=data_1.9[-5,-1]) #5번쨰 데이터 제외
#회귀계수
data.frame(all=coef(res_1),
rm4=coef(res_2),
rm5=coef(res_3)) all rm4 rm5
(Intercept) 1.722213529 1.099471134 1.626014115
Agr 0.005809126 0.010136685 0.002352222
Forest -0.012967887 -0.007589231 -0.012760349
Rsdntial -0.007226768 -0.123792917 0.181160986
ComIndl 0.305027765 1.528956204 0.075617570#p-value
data.frame(all=coef(summary(res_1))[,4],
rm4=coef(summary(res_2))[,4],
rm5=coef(summary(res_3))[,4]) all rm4 rm5
(Intercept) 0.18316946 0.2477883387 0.05619948
Agr 0.70462624 0.3717054741 0.80880737
Forest 0.36667966 0.4700975391 0.16975563
Rsdntial 0.83372002 0.0071342930 0.00112280
ComIndl 0.08230952 0.0005512227 0.51774981#4,5번째 데이터를 각각 뺴고 진행을 해보니 영향을 끼치는 값임을 알수있고
#5번째는 주거 관련해서는 부호를 바꿀 정도로 강력하다
# 단순선형회구모형
res<-lm(Nitrogen~ComIndl,data=data_1.9)그림 [4.5]
plot(Nitrogen~ComIndl,data=data_1.9,pch=19)
abline(res)
#leverage values 지레값
p_ii<-hatvalues(res)
hiegh_leverage<-ifelse(p_ii>2*2/30,data_1.9$River,"")
hiegh_leverage #높은 지레값을 가지고 있는 강의 이름을 표시(이는 보기에 편하기 위해서함) 1 2 3 4 5 6
"" "" "" "Neversink" "Hackensack" ""
7 8 9 10 11 12
"" "" "" "" "" ""
13 14 15 16 17 18
"" "" "" "" "" ""
19 20
"" "" text(data_1.9$ComIndl,data_1.9$Nitrogen-0.1,hiegh_leverage)
그림 [4.6] - 잔차의 인덱스플롯
plot(rstandard(res),pch=19) #2또는 3시그마를 넘으면 특이값이라 함 
그림 [4.6] - 지레값의 인덱스플롯
plot(p_ii,pch=19) #평균의 2배를 기준으로 비교함
abline(h=2*2/30,col="red") #이보다 높은것이 지레값이 높은것 높은 영향력을 가진것 
# 단순선형회귀모형
res<-lm(Nitrogen~ComIndl,data=data_1.9)
plot(Nitrogen~ComIndl,data=data_1.9,pch=19)
abline(res)
res<-lm(Nitrogen~ComIndl,data=data_1.9[-4,-1])
plot(Nitrogen~ComIndl,data=data_1.9[-4,-1],pch=19)
abline(res) #4번째 제외
res<-lm(Nitrogen~ComIndl,data=data_1.9[-5,-1])
plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)
abline(res) #5번쨰 제외
#이처럼 4,5번을 빼고 진행을 하면 조금더 잘 나타냄
res<-lm(Nitrogen~ComIndl,data=data_1.9[c(-4,-5),-1])
plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)
abline(res) #4,5번째 제외 
4.9 영향력의 측도
4.9.1 Cook의 거리
res<-lm(Nitrogen~ComIndl,data=data_1.9)
plot(res,4)
#install.packages("olsrr")
library(olsrr)Warning: package 'olsrr' was built under R version 4.3.2
Attaching package: 'olsrr'The following object is masked from 'package:datasets':
rivers
ols_plot_cooksd_chart(res)
4.9.2 Welsch & Kuh의 측도
olsrr::ols_plot_dffits(res)
4.9.3 Hadi의 영향력 측도
olsrr::ols_plot_hadi(res)
# Residual & Leverage & Cook's distance
plot(res,5) #영향력 관측치를 보기 위한 플랏 
4.10 잠재성 - 잔차플롯
olsrr::ols_plot_resid_pot(res) #2.0에 있는 값외에도 x축 0.2이상의 것들도 특이값으로 
#지레값이 커도 영향력이 없는 애들은 신경 안써도 되나 영향력이 큰애들을 탐색해 보아야 함4.11 특이값에 대한 처리
#특이값이 큰경우 그것이 TRUE데이터면 오류가 있는 경우 수정을 하거나 가중치를 변화를
#하거나 데이터를 수정을 시켜주거나 다시 실험을 하는 등 여러가지 방법을 사용하여..
## 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)
res<-lm(Nitrogen~ComIndl,data=data_1.9)
car::avPlots(res,pch=19)
res<-lm(Nitrogen~.,data=data_1.9[-1])
car::avPlots(res,pch=19)
#beta별로 각각의 어떤 변수가 영향력을 많이 주는지 알게하는 함수 사례:스코틀랜드 언덕 경주 데이터
data_4.5<-read.table("All_Data/p120.txt",header=T,sep="\t")
dim(data_4.5)[1] 35 4names(data_4.5)[1] "Hill.Race" "Time" "Distance" "Climb" head(data_4.5) Hill.Race Time Distance Climb
1 Greenmantle New Year Dash 965 2.5 650
2 Carnethy 2901 6.0 2500
3 Craig Dunain 2019 6.0 900
4 Ben Rha 2736 7.5 800
5 Ben Lomond 3736 8.0 3070
6 Goatfell 4393 8.0 2866# 회전도표
library(rgl)
with(data_4.5,plot3d(x=Distance,y=Climb,z=Time))
res<-lm(Time~Distance+Climb,data=data_4.5)
summary(res) #beta_0가 마이너스여도 신경안씀 관심있는 것은 회귀계수임
Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)
Residuals:
Min 1Q Median 3Q Max
-973.0 -427.7 -71.2 142.2 3907.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -539.4829 258.1607 -2.090 0.0447 *
Distance 373.0727 36.0684 10.343 9.86e-12 ***
Climb 0.6629 0.1231 5.387 6.44e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 880.5 on 32 degrees of freedom
Multiple R-squared: 0.9191, Adjusted R-squared: 0.914
F-statistic: 181.7 on 2 and 32 DF, p-value: < 2.2e-16#Time=-539.4829+373.0727*Distance+0.6629*Climb ### 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)
car::avPlots(res,pch=19)
### 성분잔차플롯(component plus residual plots), 편자차플롯(partial residual plot)
car::crPlots(res,id=T,pch=19) #첨가변수보다 성분잔차를 더 많이 사용 / 비선형여부를 확인
# 점선에 비해서 분홍선이 크게 떨어져 있지 않아 선형적인 추세를 가지고 있다고 추정이 가능
### 잠재성-잔차플롯
olsrr::ols_plot_resid_pot(res)
### Hadi의 영향력 측도
olsrr::ols_plot_hadi(res)
### Cook의 거리
olsrr::ols_plot_cooksd_chart(res) #전체적은 플롯을 보면서 이상치에 대한 확인을 함 
# 이러한 값들의 제외 여부는 연구자가 선택해서 진행을 함
# outlier에 대한 처리를 어떻게 했다고 말을 해야함
# 보고서는 옆에 사람이 보고 쉽게 따라할 수 있을 정도로
# 어떤 속성으로 어떻게 진행을 했다라는 것을 중시 코드 보단 결과를 보여줘라 4.14 로버스트 회귀 - outlier에 크게 영향을 받지 않음
library(MASS)
Attaching package: 'MASS'The following object is masked from 'package:olsrr':
cementThe following object is masked from 'package:dplyr':
selectres<-lm(Time~Distance+Climb,data=data_4.5)
res
Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)
Coefficients:
(Intercept) Distance Climb
-539.4829 373.0727 0.6629 res_rlm<-MASS::rlm(Time~Distance+Climb,data=data_4.5)
res_rlmCall:
rlm(formula = Time ~ Distance + Climb, data = data_4.5)
Converged in 10 iterations
Coefficients:
(Intercept) Distance Climb
-576.3836570 393.0374614 0.4977894
Degrees of freedom: 35 total; 32 residual
Scale estimate: 313 summary(res)
Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)
Residuals:
Min 1Q Median 3Q Max
-973.0 -427.7 -71.2 142.2 3907.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -539.4829 258.1607 -2.090 0.0447 *
Distance 373.0727 36.0684 10.343 9.86e-12 ***
Climb 0.6629 0.1231 5.387 6.44e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 880.5 on 32 degrees of freedom
Multiple R-squared: 0.9191, Adjusted R-squared: 0.914
F-statistic: 181.7 on 2 and 32 DF, p-value: < 2.2e-16summary(res_rlm)
Call: rlm(formula = Time ~ Distance + Climb, data = data_4.5)
Residuals:
Min 1Q Median 3Q Max
-645.074 -197.082 -2.035 212.266 3942.045
Coefficients:
Value Std. Error t value
(Intercept) -576.3837 105.2774 -5.4749
Distance 393.0375 14.7086 26.7216
Climb 0.4978 0.0502 9.9200
Residual standard error: 312.6 on 32 degrees of freedom#install.packages("robustbase")
library(robustbase)Warning: package 'robustbase' was built under R version 4.3.2res_lmrob<-lmrob(Time~Distance+Climb,data=data_4.5)
res_lmrob
Call:
lmrob(formula = Time ~ Distance + Climb, data = data_4.5)
\--> method = "MM"
Coefficients:
(Intercept) Distance Climb
-487.3793 398.2784 0.3901 summary(res_lmrob)
Call:
lmrob(formula = Time ~ Distance + Climb, data = data_4.5)
\--> method = "MM"
Residuals:
Min 1Q Median 3Q Max
-650.01 -160.60 23.37 216.11 3875.00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -487.37929 86.33384 -5.645 3.04e-06 ***
Distance 398.27843 5.98522 66.544 < 2e-16 ***
Climb 0.39013 0.04165 9.368 1.09e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Robust residual standard error: 290.8
Multiple R-squared: 0.9868, Adjusted R-squared: 0.986
Convergence in 9 IRWLS iterations
Robustness weights:
3 observations c(7,18,33) are outliers with |weight| = 0 ( < 0.0029);
2 weights are ~= 1. The remaining 30 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.5965 0.9098 0.9623 0.9214 0.9875 0.9990
Algorithmic parameters:
tuning.chi bb tuning.psi refine.tol
1.548e+00 5.000e-01 4.685e+00 1.000e-07
rel.tol scale.tol solve.tol zero.tol
1.000e-07 1.000e-10 1.000e-07 1.000e-10
eps.outlier eps.x warn.limit.reject warn.limit.meanrw
2.857e-03 1.364e-08 5.000e-01 5.000e-01
nResample max.it best.r.s k.fast.s k.max
500 50 2 1 200
maxit.scale trace.lev mts compute.rd fast.s.large.n
200 0 1000 0 2000
psi subsampling cov
"bisquare" "nonsingular" ".vcov.avar1"
compute.outlier.stats
"SM"
seed : int(0) #standard error가 res > res_rlm > res_lmrob 순으로 되어 있음
# 4장 연습문제 해보기0509 - regression analysis
library(dplyr)5장 질적 예측 변수
#install.packages("fastDummies")
library(fastDummies)Warning: package 'fastDummies' was built under R version 4.3.2Thank you for using fastDummies!To acknowledge our work, please cite the package:Kaplan, J. & Schlegel, B. (2023). fastDummies: Fast Creation of Dummy (Binary) Columns and Rows from Categorical Variables. Version 1.7.1. URL: https://github.com/jacobkap/fastDummies, https://jacobkap.github.io/fastDummies/.5.2 급료조사 데이터
data_5.1<-read.table("All_Data/p130.txt",header=T,sep="\t")
# S:급료 X:경력 E:교육수준 M:관리(형태) / E,M은 범주형 변수
names(data_5.1)[1] "S" "X" "E" "M"# 범주형 질적 변수를 수치형으로 변형시켜서 예측하는데 사용한것이 질적 예측변수이다
# E_1,E_2,E_3이런식으로 나누어서 0,1로 분류를 한다 (이것이 가변수)
# 더미변수를 만들 경우에는 역행렬의 조건에 의해서 -1개의 변수만 만들면 된다
# 이는 공산성의 문제또한 있기에 이를 위해서 -1를 한것임
### 자료형 변경 : 정수 -> 범주
data_5.1$E<-as.factor(data_5.1$E)
data_5.1$M<-as.factor(data_5.1$M)
head(data_5.1) S X E M
1 13876 1 1 1
2 11608 1 3 0
3 18701 1 3 1
4 11283 1 2 0
5 11767 1 3 0
6 20872 2 2 1data_5.1$E #Levels: 1 2 3이라는 것이 생김 문자로 처리한다는 의미 [1] 1 3 3 2 3 2 2 1 3 2 1 2 3 1 3 3 2 2 3 1 1 3 2 2 1 2 1 3 1 1 2 3 2 2 1 2 3 1
[39] 2 2 3 2 2 1 2 1
Levels: 1 2 3### 가변수 생성
data_5.1$E<-factor(as.character(data_5.1$E),levels = c("3","1","2"))
#3번을 베이스 카테고리로 쓰기위한 설정 / 설정을 안하면 베이스는 e_1이 된다
data_5.1$M<-factor(as.character(data_5.1$M),levels = c("0","1"))
data_dummy<-dummy_cols(data_5.1,
select_columns = c("E","M"),
remove_first_dummy = T,
remove_selected_columns = T) #첫번째 생성되는 더미 변수를 제거
data_dummy #가변수의 더미는 n-1개를 하는 것이 역행렬을 위한 것이기에 지워준다 S X E_1 E_2 M_1
1 13876 1 1 0 1
2 11608 1 0 0 0
3 18701 1 0 0 1
4 11283 1 0 1 0
5 11767 1 0 0 0
6 20872 2 0 1 1
7 11772 2 0 1 0
8 10535 2 1 0 0
9 12195 2 0 0 0
10 12313 3 0 1 0
11 14975 3 1 0 1
12 21371 3 0 1 1
13 19800 3 0 0 1
14 11417 4 1 0 0
15 20263 4 0 0 1
16 13231 4 0 0 0
17 12884 4 0 1 0
18 13245 5 0 1 0
19 13677 5 0 0 0
20 15965 5 1 0 1
21 12336 6 1 0 0
22 21352 6 0 0 1
23 13839 6 0 1 0
24 22884 6 0 1 1
25 16978 7 1 0 1
26 14803 8 0 1 0
27 17404 8 1 0 1
28 22184 8 0 0 1
29 13548 8 1 0 0
30 14467 10 1 0 0
31 15942 10 0 1 0
32 23174 10 0 0 1
33 23780 10 0 1 1
34 25410 11 0 1 1
35 14861 11 1 0 0
36 16882 12 0 1 0
37 24170 12 0 0 1
38 15990 13 1 0 0
39 26330 13 0 1 1
40 17949 14 0 1 0
41 25685 15 0 0 1
42 27837 16 0 1 1
43 18838 16 0 1 0
44 17483 16 1 0 0
45 19207 17 0 1 0
46 19346 20 1 0 0# 더미를 만든 이후 더미의 모체 변수인 E M을 지워주어야 한다
### 회귀분석(1) - 가변수
res<-lm(S~.,data = data_dummy)
res
Call:
lm(formula = S ~ ., data = data_dummy)
Coefficients:
(Intercept) X E_1 E_2 M_1
11031.8 546.2 -2996.2 147.8 6883.5 ### 회귀분석(2) - lm()
res_1<-lm(S~.,data=data_5.1)
res_1
Call:
lm(formula = S ~ ., data = data_5.1)
Coefficients:
(Intercept) X E1 E2 M1
11031.8 546.2 -2996.2 147.8 6883.5 #더미변수를 많이 쓰기에 factor로 바꾸어주고 분석하면 알아서 더미변수를 만들어서 진행함
summary(res_1)
Call:
lm(formula = S ~ ., data = data_5.1)
Residuals:
Min 1Q Median 3Q Max
-1884.60 -653.60 22.23 844.85 1716.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11031.81 383.22 28.787 < 2e-16 ***
X 546.18 30.52 17.896 < 2e-16 ***
E1 -2996.21 411.75 -7.277 6.72e-09 ***
E2 147.82 387.66 0.381 0.705
M1 6883.53 313.92 21.928 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1027 on 41 degrees of freedom
Multiple R-squared: 0.9568, Adjusted R-squared: 0.9525
F-statistic: 226.8 on 4 and 41 DF, p-value: < 2.2e-16#E_3와 M_0는 Intercept에 포함이 되어있음 그렇기에 베이스 카테고리라고 한다
#E가 3이고 M이 0이면 대학원이상 일반관리 직급 -> Intercept(Beta_0) + X*Beta_1 = Y
### 표준화잔차 대 경력연수
plot(data_5.1$X, rstandard(res_1),pch=19,xlab="범주",ylab="잔차")
# 0을 중심으로 잘 퍼져있는가를 확인해야함
#### 표준화잔차 대 교육수준 - 관리 조합
EM<-paste0(data_5.1$E,data_5.1$M)
plot(EM,rstandard(res_1),pch=19,xlabl="범주",ylab="잔차")Warning in plot.window(...): "xlabl" is not a graphical parameterWarning in plot.xy(xy, type, ...): "xlabl" is not a graphical parameterWarning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameterWarning in box(...): "xlabl" is not a graphical parameterWarning in title(...): "xlabl" is not a graphical parameter
### 상호작용 효과(Interaction Effect)
res<-lm(S~X+E+M+E*M,data=data_5.1)
res
Call:
lm(formula = S ~ X + E + M + E * M, data = data_5.1)
Coefficients:
(Intercept) X E1 E2 M1 E1:M1
11203.4 497.0 -1730.7 -349.1 7047.4 -3066.0
E2:M1
1836.5 summary(res)
Call:
lm(formula = S ~ X + E + M + E * M, data = data_5.1)
Residuals:
Min 1Q Median 3Q Max
-928.13 -46.21 24.33 65.88 204.89
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11203.434 79.065 141.698 < 2e-16 ***
X 496.987 5.566 89.283 < 2e-16 ***
E1 -1730.748 105.334 -16.431 < 2e-16 ***
E2 -349.078 97.568 -3.578 0.000945 ***
M1 7047.412 102.589 68.695 < 2e-16 ***
E1:M1 -3066.035 149.330 -20.532 < 2e-16 ***
E2:M1 1836.488 131.167 14.001 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 173.8 on 39 degrees of freedom
Multiple R-squared: 0.9988, Adjusted R-squared: 0.9986
F-statistic: 5517 on 6 and 39 DF, p-value: < 2.2e-16### 표준화잔차 대 경력연수
plot(data_5.1$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")
### Cook의 거리
plot(res,4) #33번째 데이터만 제외하고 다시 회귀모형을 생성예정
### 상호작용 효과 - 관측개체 33 제외
data_use<-data_5.1[-33,]
res<-lm(S~X+E+M+E*M,data=data_use)
res
Call:
lm(formula = S ~ X + E + M + E * M, data = data_use)
Coefficients:
(Intercept) X E1 E2 M1 E1:M1
11199.7 498.4 -1741.3 -357.0 7040.6 -3051.8
E2:M1
1997.5 summary(res)
Call:
lm(formula = S ~ X + E + M + E * M, data = data_use)
Residuals:
Min 1Q Median 3Q Max
-112.884 -43.636 -5.036 46.622 128.480
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11199.714 30.533 366.802 < 2e-16 ***
X 498.418 2.152 231.640 < 2e-16 ***
E1 -1741.336 40.683 -42.803 < 2e-16 ***
E2 -357.042 37.681 -9.475 1.49e-11 ***
M1 7040.580 39.619 177.707 < 2e-16 ***
E1:M1 -3051.763 57.674 -52.914 < 2e-16 ***
E2:M1 1997.531 51.785 38.574 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 67.12 on 38 degrees of freedom
Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
F-statistic: 3.543e+04 on 6 and 38 DF, p-value: < 2.2e-16### 표준화잔차 대 경력연수
plot(data_use$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")
#### 표준화잔차 대 교육수준 - 관리 조합
EM<-paste0(data_use$E,data_use$M)
plot(EM,rstandard(res),pch=19,xlabl="범주",ylab="잔차")Warning in plot.window(...): "xlabl" is not a graphical parameterWarning in plot.xy(xy, type, ...): "xlabl" is not a graphical parameterWarning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameterWarning in box(...): "xlabl" is not a graphical parameterWarning in title(...): "xlabl" is not a graphical parameter
#상호 호과가 들어간 이 모형이 더 괜찮은 모양이라고 판단이 된다
### 기본 급료 추정 - 표 5.6
res<-lm(S~X+E+M+E*M,data=data_use)
df_new<-data.frame(X=rep(0,6),
E=rep(1:3,c(2,2,2)),
M=rep(c(0,1),3))
### 가변수 생성 - 분석용
df_new$E<-factor(as.character(df_new$E),levels = c("3","1","2"))
df_new$M<-factor(as.character(df_new$M),levels = c("0","1"))
cbind(df_new,predict = predict(res,df_new,interval = "confidence")) X E M predict.fit predict.lwr predict.upr
1 0 1 0 9458.378 9395.539 9521.216
2 0 1 1 13447.195 13382.933 13511.456
3 0 2 0 10842.672 10789.719 10895.624
4 0 2 1 19880.782 19814.090 19947.474
5 0 3 0 11199.714 11137.902 11261.525
6 0 3 1 18240.294 18182.503 18298.084# 이를 보고 각 레벨에 따른 차이를 보고 얼마나 나는 지 분석이 가능해야 한다
# ex) 고졸과 대학원졸의 관리자 직급의 급여의 차이는?(평균적으로)5.4 회귀방정식의 체계: 두 집단의 비교
표 5.7 - 고용전 검사 프로그램 데이터
data_5.7<-read.table("All_Data/p140.txt",header=T,sep="\t")
head(data_5.7) TEST RACE JPERF
1 0.28 1 1.83
2 0.97 1 4.59
3 1.25 1 2.97
4 2.46 1 8.14
5 2.51 1 8.00
6 1.17 1 3.30# 모형 1 - 통합모형 인종간 차이가 없을때
model_1<-lm(JPERF~TEST,data=data_5.7)
summary(model_1)
Call:
lm(formula = JPERF ~ TEST, data = data_5.7)
Residuals:
Min 1Q Median 3Q Max
-3.3558 -0.8798 -0.1897 1.2735 2.3312
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0350 0.8680 1.192 0.248617
TEST 2.3605 0.5381 4.387 0.000356 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.591 on 18 degrees of freedom
Multiple R-squared: 0.5167, Adjusted R-squared: 0.4899
F-statistic: 19.25 on 1 and 18 DF, p-value: 0.0003555# 모형 3
model_3<-lm(JPERF~TEST+RACE+TEST*RACE,data=data_5.7)
summary(model_3)
Call:
lm(formula = JPERF ~ TEST + RACE + TEST * RACE, data = data_5.7)
Residuals:
Min 1Q Median 3Q Max
-2.0734 -1.0594 -0.2548 1.2830 2.1980
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0103 1.0501 1.914 0.0736 .
TEST 1.3134 0.6704 1.959 0.0677 .
RACE -1.9132 1.5403 -1.242 0.2321
TEST:RACE 1.9975 0.9544 2.093 0.0527 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.407 on 16 degrees of freedom
Multiple R-squared: 0.6643, Adjusted R-squared: 0.6013
F-statistic: 10.55 on 3 and 16 DF, p-value: 0.0004511# 인종적인 차이가 있는지 없는지를 확인해야 한다 어떤 모형을 사용할지
### H_0:gamma=delta=0
anova(model_1,model_3) #model_3가 FM(완전모형)Analysis of Variance Table
Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + RACE + TEST * RACE
Res.Df RSS Df Sum of Sq F Pr(>F)
1 18 45.568
2 16 31.655 2 13.913 3.5161 0.05424 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#P-value(0.05424)<0.05이기에 H_0 이기에 모형1을 선택하는 것이 옳다고 판단(그러나 확신x)
#서로의 R-squared롤 보면 model_3가 더 좋음 / ANOVA는 참고용 절대적이지는 않다 그림 5.7 - 표준화잔차 대 검사점수 : 모형 1
plot(data_5.7$TEST,rstandard(model_1),
pch=19,xlab="검사점수",ylab="잔차",
main="그림 5.7 - 표준화잔차 대 검사점수 : 모형 1")
그림 5.8 - 표준화잔차 대 검사점수 : 모형 3
plot(data_5.7$TEST,rstandard(model_3),
pch=19,xlab="검사점수",ylab="잔차",
main="그림 5.8 - 표준화잔차 대 검사점수 : 모형 3")
# 통계에서 나오는 결과는 결정의 보조수단이지 절대적이지 않아
# 3번 모형을 선택한다고 결정한다고 진행
summary(model_3) # Multiple R-squared: 0.6643
Call:
lm(formula = JPERF ~ TEST + RACE + TEST * RACE, data = data_5.7)
Residuals:
Min 1Q Median 3Q Max
-2.0734 -1.0594 -0.2548 1.2830 2.1980
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0103 1.0501 1.914 0.0736 .
TEST 1.3134 0.6704 1.959 0.0677 .
RACE -1.9132 1.5403 -1.242 0.2321
TEST:RACE 1.9975 0.9544 2.093 0.0527 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.407 on 16 degrees of freedom
Multiple R-squared: 0.6643, Adjusted R-squared: 0.6013
F-statistic: 10.55 on 3 and 16 DF, p-value: 0.0004511그림 5.9 - 표준화잔차 대 검사점수 : 모형 1
plot(data_5.7$RACE,rstandard(model_1),
pch=19,xlab="검사점수",ylab="잔차",
main="그림 5.9 - 표준화잔차 대 검사점수 : 모형 1")
# 분리된 회귀분석 결과
data_5.7_R1<-subset(data_5.7,RACE==1)
model_R1<-lm(JPERF~TEST,data=data_5.7_R1)
summary(model_R1) #소수민족
Call:
lm(formula = JPERF ~ TEST, data = data_5.7_R1)
Residuals:
Min 1Q Median 3Q Max
-2.0734 -0.6267 -0.2548 1.1624 1.5394
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09712 1.03519 0.094 0.927564
TEST 3.31095 0.62411 5.305 0.000724 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.292 on 8 degrees of freedom
Multiple R-squared: 0.7787, Adjusted R-squared: 0.751
F-statistic: 28.14 on 1 and 8 DF, p-value: 0.0007239data_5.7_R0<-subset(data_5.7,RACE==0)
model_R0<-lm(JPERF~TEST,data=data_5.7_R0)
summary(model_R0) #백인
Call:
lm(formula = JPERF ~ TEST, data = data_5.7_R0)
Residuals:
Min 1Q Median 3Q Max
-1.8599 -1.0663 -0.3061 1.0957 2.1980
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0103 1.1291 1.780 0.113
TEST 1.3134 0.7208 1.822 0.106
Residual standard error: 1.512 on 8 degrees of freedom
Multiple R-squared: 0.2933, Adjusted R-squared: 0.205
F-statistic: 3.32 on 1 and 8 DF, p-value: 0.1059# 통합모형에서 나온 각각의 회귀식이 통합모형에서 나온것과 동일함 따라서 인종별로 나누어서
# 진행할 필요없이 통일모형을 사용해서 진행을 하면 된다(이는 데이터셋을 나눈경우와 동일함)그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만
plot(data_5.7_R1$TEST,rstandard(model_R1),
pch=19,xlab="검사점수",ylab="잔차",
main="그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만만")
그림 5.11 - 표준화잔차 대 검사점수 : 모형 1 . 백인만
plot(data_5.7_R0$TEST,rstandard(model_R0),
pch=19,xlab="검사점수",ylab="잔차",
main="그림 5.11 - 표준화잔차 대 검사점수 : 모형 1. 백인만")
# 적절한 합격점수의 결정 - 소수민족
# 고용전 검사점수의 합격점에 대한 95% 신뢰구간
ym<-4
xm<-(ym-0.09712)/3.31095
s<-1.292
n<-10
t<-qt(1-0.05/2,8)
c(xm-(t*s/n)/3.31095, xm+(t*s/n)/3.31095) #신뢰구간[1] 1.088795 1.2687645.4.2 동일한 기울기와 다른 절편항을 가지는 모형
model_3<-lm(JPERF~TEST+RACE,data=data_5.7)
summary(model_3)
Call:
lm(formula = JPERF ~ TEST + RACE, data = data_5.7)
Residuals:
Min 1Q Median 3Q Max
-2.7872 -1.0370 -0.2095 0.9198 2.3645
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6120 0.8870 0.690 0.499578
TEST 2.2988 0.5225 4.400 0.000391 ***
RACE 1.0276 0.6909 1.487 0.155246
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.54 on 17 degrees of freedom
Multiple R-squared: 0.5724, Adjusted R-squared: 0.5221
F-statistic: 11.38 on 2 and 17 DF, p-value: 0.0007312#Intercept = BETA_0 / TEST = BETA_1 / RACE = gamma
## H_0: gamma=0
anova(model_1,model_3) #p-value(0.1552)>0.05 이기에 H_0은 참이다//gamma=0Analysis of Variance Table
Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + RACE
Res.Df RSS Df Sum of Sq F Pr(>F)
1 18 45.568
2 17 40.322 1 5.2468 2.2121 0.1552#기울기가 같고 절편이 다른 모형은 아니라고 데이터가 이야기하고 있다
# 소수민족(RACE=1): (0.6120+1.0276)+2.2988*TEST = 1.6396+2.2988*TEST
# 백인(RACE=0): (0.6120)+2.2988*TEST5.4.3 동일한 절편항과 다른 기울기를 가지는 모형
model_3<-lm(JPERF~TEST+I(TEST*RACE),data=data_5.7)
summary(model_3) #I()를 하면 교호작용하는 것만 보이게 하려고 없으면 RACE항이 자동추가됨
Call:
lm(formula = JPERF ~ TEST + I(TEST * RACE), data = data_5.7)
Residuals:
Min 1Q Median 3Q Max
-2.41100 -0.88871 -0.03359 0.97720 2.44440
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1211 0.7804 1.437 0.16900
TEST 1.8276 0.5356 3.412 0.00332 **
I(TEST * RACE) 0.9161 0.3972 2.306 0.03395 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.429 on 17 degrees of freedom
Multiple R-squared: 0.6319, Adjusted R-squared: 0.5886
F-statistic: 14.59 on 2 and 17 DF, p-value: 0.0002045## H_0: delta=0
anova(model_1,model_3) #p-value(0.03395)<0.05보다 작기에 delta항은 필요한 변수라는 사실 Analysis of Variance Table
Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + I(TEST * RACE)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 18 45.568
2 17 34.708 1 10.861 5.3196 0.03395 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Intercept = BETA_0 / TEST = BETA_1 / I(TEST*RACE) = delta6장 변수변환 : Transformation of Varibables
library(dplyr)6.3 x-선 방사에 의한 박테리아 사망률
data_6.2<-read.table("All_Data/p168.txt",header=T,sep="\t")
head(data_6.2) t N_t
1 1 355
2 2 211
3 3 197
4 4 166
5 5 142
6 6 106그림 6.5
plot(N_t~t,data=data_6.2,pch=19)
6.3.1 선형모형의 부적절성
res<-lm(N_t~t,data=data_6.2)
summary(res)
Call:
lm(formula = N_t ~ t, data = data_6.2)
Residuals:
Min 1Q Median 3Q Max
-43.867 -23.599 -9.652 10.223 114.883
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 259.58 22.73 11.420 3.78e-08 ***
t -19.46 2.50 -7.786 3.01e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 41.83 on 13 degrees of freedom
Multiple R-squared: 0.8234, Adjusted R-squared: 0.8098
F-statistic: 60.62 on 1 and 13 DF, p-value: 3.006e-06그림 6.6
plot(data_6.2$t,rstandard(res),pch=19) # 표준화 잔차의 플랏
# 잔차가 랜덤하게 잘 퍼져있어야 하는데 적절한 회귀모형이 아니라는 사실이 나옴6.3.2 선형성을 위한 로그변환
그림 6.7
plot(log(N_t)~t,data=data_6.2,pch=19)
# 박테리아의 수에 로그를 취하니 선형성이 보인다
res<-lm(log(N_t)~t,data=data_6.2)
summary(res)
Call:
lm(formula = log(N_t) ~ t, data = data_6.2)
Residuals:
Min 1Q Median 3Q Max
-0.18445 -0.06189 0.01253 0.05201 0.20021
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.973160 0.059778 99.92 < 2e-16 ***
t -0.218425 0.006575 -33.22 5.86e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.11 on 13 degrees of freedom
Multiple R-squared: 0.9884, Adjusted R-squared: 0.9875
F-statistic: 1104 on 1 and 13 DF, p-value: 5.86e-14plot(data_6.2$t,rstandard(res),pch=19)
# 로그를 취해주니 잔차가 랜덤하게 이루어져 있음
# n_0에 대한 추론
exp(5.973160)[1] 392.7448exp(coef(res)[1]) #로그를 취해주고 하는 부분이 잘 이해가 안됨 (Intercept)
392.7449 exp(coef(res)[1]-0.0588/2) #불편추정량 구하기 (참고용)(Intercept)
381.3663 6.4 분산안정화 변환
data_6.6<-read.table("All_Data/p174.txt",header=T,sep="\t")
data_6.6 #운항률과 사고 발생 건수에 대한 자료 Y N
1 11 0.0950
2 7 0.1920
3 7 0.0750
4 19 0.2078
5 9 0.1382
6 4 0.0540
7 3 0.1292
8 1 0.0503
9 3 0.0629plot(Y~N,data=data_6.6,pch=19) #잔차의 분산이 계속 커지는 효과를 보임 
res_1<-lm(Y~N,data=data_6.6)
summary(res_1)
Call:
lm(formula = Y ~ N, data = data_6.6)
Residuals:
Min 1Q Median 3Q Max
-5.3351 -2.1281 0.1605 2.2670 5.6382
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1402 3.1412 -0.045 0.9657
N 64.9755 25.1959 2.579 0.0365 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.201 on 7 degrees of freedom
Multiple R-squared: 0.4872, Adjusted R-squared: 0.4139
F-statistic: 6.65 on 1 and 7 DF, p-value: 0.03654# 표준화잔차 대 N의 플롯 / 그림 6.11
plot(data_6.6$N,rstandard(res_1),pch=19) #등분산성에 만족하지 못하는 모형을 보인다 
# N에 대한 sqrt(Y)의 회귀
# N에 대한 Y의 회귀
res_2<-lm(sqrt(Y)~N,data=data_6.6)
summary(res_2)
Call:
lm(formula = sqrt(Y) ~ N, data = data_6.6)
Residuals:
Min 1Q Median 3Q Max
-0.9690 -0.7655 0.1906 0.5874 1.0211
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1692 0.5783 2.022 0.0829 .
N 11.8564 4.6382 2.556 0.0378 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7733 on 7 degrees of freedom
Multiple R-squared: 0.4828, Adjusted R-squared: 0.4089
F-statistic: 6.535 on 1 and 7 DF, p-value: 0.03776그림 6.12
plot(data_6.6$N,rstandard(res_2),pch=19) #0을 중심으로 잔차가 보다 잘 퍼져있음이 보임 
6.5 이분산성의 검출
data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")
data_6.9 X Y
1 294 30
2 247 32
3 267 37
4 358 44
5 423 47
6 311 49
7 450 56
8 534 62
9 438 68
10 697 78
11 688 80
12 630 84
13 709 88
14 627 97
15 615 100
16 999 109
17 1022 114
18 1015 117
19 700 106
20 850 128
21 980 130
22 1025 160
23 1021 97
24 1200 180
25 1250 112
26 1500 210
27 1650 135# Y 대 X의 플로
plot(Y~X,data=data_6.9,pch=19,main="그림 6.13")
# X에 대한 Y의 회귀
res_1<-lm(Y~X,data=data_6.9)
summary(res_1)
Call:
lm(formula = Y ~ X, data = data_6.9)
Residuals:
Min 1Q Median 3Q Max
-53.294 -9.298 -5.579 14.394 39.119
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.44806 9.56201 1.511 0.143
X 0.10536 0.01133 9.303 1.35e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.73 on 25 degrees of freedom
Multiple R-squared: 0.7759, Adjusted R-squared: 0.7669
F-statistic: 86.54 on 1 and 25 DF, p-value: 1.35e-09# 표준화잔차 대 X의 플롯
plot(data_6.9$X, rstandard(res_1),pch=19,main = "그림 6.14")
6.6 이분산성의 제거
# 변환된 Y/X와 1/X를 적합한 회귀
data_6.9_1<-data.frame(Y=data_6.9$Y/data_6.9$X,
X=1/data_6.9$X)
res_2<-lm(Y~X,data=data_6.9_1)
summary(res_2) #지금은 변환하고 프라임 값들의 추정치가 나온것이기에 원래 회귀식으로 돌아가야함
Call:
lm(formula = Y ~ X, data = data_6.9_1)
Residuals:
Min 1Q Median 3Q Max
-0.041477 -0.013852 -0.004998 0.024671 0.035427
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.120990 0.008999 13.445 6.04e-13 ***
X 3.803296 4.569745 0.832 0.413
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.02266 on 25 degrees of freedom
Multiple R-squared: 0.02696, Adjusted R-squared: -0.01196
F-statistic: 0.6927 on 1 and 25 DF, p-value: 0.4131# 본래 변환시 X를 나눴기에 X를 곱해준다 -> B_0,B_1의 추정값이 바뀐다 서로
plot(data_6.9_1$X, rstandard(res_2),pch=19,
xlab="1/X",ylab="잔차",main = "[그림 6.15]")
6.7 가중최소제곱법
wt<-1/data_6.9$X^2 #가중값
res_3<-lm(Y~X,data=data_6.9,weights = wt)
summary(res_3) #결과는 위의 6.6과 동일한 값이 나옴
Call:
lm(formula = Y ~ X, data = data_6.9, weights = wt)
Weighted Residuals:
Min 1Q Median 3Q Max
-0.041477 -0.013852 -0.004998 0.024671 0.035427
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.803296 4.569745 0.832 0.413
X 0.120990 0.008999 13.445 6.04e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.02266 on 25 degrees of freedom
Multiple R-squared: 0.8785, Adjusted R-squared: 0.8737
F-statistic: 180.8 on 1 and 25 DF, p-value: 6.044e-136.8 데이터에 대한 로그변환
data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")
head(data_6.9) X Y
1 294 30
2 247 32
3 267 37
4 358 44
5 423 47
6 311 49# log(Y) 대 X의 산점도
plot((Y)~X,data=data_6.9,pch=19,main="식 6.9")
plot(log(Y)~X,data=data_6.9,pch=19,main="그림 6.16")
res_4<-lm(log(Y)~X,data=data_6.9)
summary(res_4)
Call:
lm(formula = log(Y) ~ X, data = data_6.9)
Residuals:
Min 1Q Median 3Q Max
-0.59648 -0.16578 0.00244 0.17481 0.34964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.5150232 0.1110670 31.648 < 2e-16 ***
X 0.0012041 0.0001316 9.153 1.85e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2524 on 25 degrees of freedom
Multiple R-squared: 0.7702, Adjusted R-squared: 0.761
F-statistic: 83.77 on 1 and 25 DF, p-value: 1.855e-09# X에 대한 log(Y)의 회귀로 부터 얻은 표준화잔차플롯
plot(data_6.9$X,rstandard(res_4),pch=19,
xlab="X",ylab="잔차",main="그림 6.17")
# X와 X^2에 대한 log(Y)의 회귀
df<-data.frame(log_Y=log(data_6.9$Y),
X=data_6.9$X,
X2=data_6.9$X^2)
res_5<-lm(log_Y~X+X2,data=df) #그러나 이러한 과정은 필요업고 아래방식으로..
res_5<-lm(log(Y)~X+I(X^2),data=data_6.9)
summary(res_5)
Call:
lm(formula = log(Y) ~ X + I(X^2), data = data_6.9)
Residuals:
Min 1Q Median 3Q Max
-0.30589 -0.11705 -0.02707 0.17593 0.30657
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.852e+00 1.566e-01 18.205 1.50e-15 ***
X 3.113e-03 3.989e-04 7.803 4.90e-08 ***
I(X^2) -1.102e-06 2.238e-07 -4.925 5.03e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1817 on 24 degrees of freedom
Multiple R-squared: 0.8857, Adjusted R-squared: 0.8762
F-statistic: 92.98 on 2 and 24 DF, p-value: 4.976e-12# 표준화잔차 대 적합값 플롯 : X와 X^2에 대한 log(Y)의 회귀
plot(fitted(res_5),rstandard(res_5),pch=19,
xlab="X",ylab="잔차",main="그림 6.18")
# 표준화잔차 대 X의 플롯 : X와 X^2에 대한 log(Y)의 회귀
plot(data_6.9$X,rstandard(res_5),pch=19,
xlab="X",ylab="잔차",main="그림 6.19")
# 표준화잔차 대 X^2의 플롯 : X와 X^2에 대한 log(Y)의 회귀
plot(data_6.9$X^2,rstandard(res_5),pch=19,
xlab="X",ylab="잔차",main="그림 6.20") #2차항을 추가함으로서 더 잘 피팅됨 
# 7장은 스킵 / 8장 스킵9.2 통계적 추론에 미치는 효과
data_9.1<-read.table("All_Data/p236.txt",header=T,sep="\t")
head(data_9.1) ACHV FAM PEER SCHOOL
1 -0.43148 0.60814 0.03509 0.16607
2 0.79969 0.79369 0.47924 0.53356
3 -0.92467 -0.82630 -0.61951 -0.78635
4 -2.19081 -1.25310 -1.21675 -1.04076
5 -2.84818 0.17399 -0.18517 0.14229
6 -0.66233 0.20246 0.12764 0.27311res<-lm(ACHV~.,data_9.1)
summary(res)
Call:
lm(formula = ACHV ~ ., data = data_9.1)
Residuals:
Min 1Q Median 3Q Max
-5.2096 -1.3934 -0.2947 1.1415 4.5881
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06996 0.25064 -0.279 0.781
FAM 1.10126 1.41056 0.781 0.438
PEER 2.32206 1.48129 1.568 0.122
SCHOOL -2.28100 2.22045 -1.027 0.308
Residual standard error: 2.07 on 66 degrees of freedom
Multiple R-squared: 0.2063, Adjusted R-squared: 0.1702
F-statistic: 5.717 on 3 and 66 DF, p-value: 0.001535#F-statistic: 5.717 on 3 and 66 DF, p-value: 0.001535
#이는 의미있는 beta가 존재한다라는 의미
#그러나 각 계수들의 회귀계수를 보니 모두 0이라는 결과가 나옴 이는 다중공선성이다
# Correlation panel
panel.cor<-function(x,y){
par(usr=c(0,1,0,1))
r<-round(cor(x,y),digits = 3)
text(0.5,0.5,r,cex=1.5)
}
pairs(data_9.1[-1],lower.panel = panel.cor)
그림 9.1
plot(fitted(res),rstandard(res),pch=19,
xlab="예측값",ylab="잔차",main="[그림 9.1]")
#이는 회귀모형에 동시에 들어가면 안된다라는 의미 (제거가 필요)9.3 예측에 미치는 효과
data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")
head(data_9.5) YEAR IMPORT DOPROD STOCK CONSUM
1 49 15.9 149.3 4.2 108.1
2 50 16.4 161.2 4.1 114.8
3 51 19.0 171.5 3.1 123.2
4 52 19.1 175.5 3.1 126.9
5 53 18.8 180.8 1.1 132.1
6 54 20.4 190.7 2.2 137.7### 산점도
pairs(data_9.5[-1],lower.panel = panel.cor)
# 회귀분석(1) : 데이터 1949~1966
res<-lm(IMPORT~.,data_9.5[-1])
summary(res) #다중공선성이 존재한다고 예측을 일단함
Call:
lm(formula = IMPORT ~ ., data = data_9.5[-1])
Residuals:
Min 1Q Median 3Q Max
-2.7208 -1.8354 -0.3479 1.2973 4.1008
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -19.7251 4.1253 -4.782 0.000293 ***
DOPROD 0.0322 0.1869 0.172 0.865650
STOCK 0.4142 0.3223 1.285 0.219545
CONSUM 0.2427 0.2854 0.851 0.409268
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.258 on 14 degrees of freedom
Multiple R-squared: 0.973, Adjusted R-squared: 0.9673
F-statistic: 168.4 on 3 and 14 DF, p-value: 3.212e-11그림 9.3 : 표준화잔차의 인덱스플롯
plot(1:nrow(data_9.5),rstandard(res),
pch=19,type="b",
xla="번호",ylab="잔차",main="[그림 9.3]")
# 회귀분석(2) : 데이터 1949~1959
data_use<-subset(data_9.5,YEAR<=59)
res<-lm(IMPORT~.,data_use[-1])
summary(res)
Call:
lm(formula = IMPORT ~ ., data = data_use[-1])
Residuals:
Min 1Q Median 3Q Max
-0.52367 -0.38953 0.05424 0.22644 0.78313
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.12799 1.21216 -8.355 6.9e-05 ***
DOPROD -0.05140 0.07028 -0.731 0.488344
STOCK 0.58695 0.09462 6.203 0.000444 ***
CONSUM 0.28685 0.10221 2.807 0.026277 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4889 on 7 degrees of freedom
Multiple R-squared: 0.9919, Adjusted R-squared: 0.9884
F-statistic: 285.6 on 3 and 7 DF, p-value: 1.112e-07그림 9.4 : 잔차플롯
plot(1:nrow(data_use),rstandard(res),
pch=19,type="b",
xla="번호",ylab="잔차",main="[그림 9.4]")
9.4 다중공선성의 탐색
# 분산확대인자
library(olsrr)
# 교육기회 균등(EEO)
res_9.1<-lm(ACHV~.,data_9.1)
summary(res_9.1)
Call:
lm(formula = ACHV ~ ., data = data_9.1)
Residuals:
Min 1Q Median 3Q Max
-5.2096 -1.3934 -0.2947 1.1415 4.5881
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06996 0.25064 -0.279 0.781
FAM 1.10126 1.41056 0.781 0.438
PEER 2.32206 1.48129 1.568 0.122
SCHOOL -2.28100 2.22045 -1.027 0.308
Residual standard error: 2.07 on 66 degrees of freedom
Multiple R-squared: 0.2063, Adjusted R-squared: 0.1702
F-statistic: 5.717 on 3 and 66 DF, p-value: 0.001535#car::vif(res_9.1)
olsrr::ols_vif_tol(res_9.1) # 10보다 크면 유의한 영향을 준다 제거필요? Variables Tolerance VIF
1 FAM 0.02660945 37.58064
2 PEER 0.03309981 30.21166
3 SCHOOL 0.01202567 83.15544표 9.5 : 프랑스 경제 데이터
data_use<-subset(data_9.5,YEAR<=59)
res_9.5<-lm(IMPORT~.,data_use[-1])
#car::vif(res_9.1)
olsrr::ols_vif_tol(res_9.5) #작으면 새로운 변수를 만들거나 제거를 통해서 진행 Variables Tolerance VIF
1 DOPROD 0.005376417 185.997470
2 STOCK 0.981441657 1.018909
3 CONSUM 0.005373166 186.110015# 상태 지수
cnd.idx<-olsrr::ols_eigen_cindex(res_9.1)
round(cnd.idx,4) Eigenvalue Condition Index intercept FAM PEER SCHOOL
1 2.9547 1.0000 0.0005 0.0030 0.0037 0.0014
2 0.9974 1.7211 0.9756 0.0000 0.0000 0.0000
3 0.0400 8.5996 0.0004 0.3068 0.4428 0.0008
4 0.0079 19.2826 0.0235 0.6903 0.5535 0.9978cnd.idx<-olsrr::ols_eigen_cindex(res_9.5)
round(cnd.idx,4) Eigenvalue Condition Index intercept DOPROD STOCK CONSUM
1 3.8384 1.0000 0.0010 0.0000 0.0109 0.0000
2 0.1484 5.0863 0.0053 0.0001 0.9385 0.0001
3 0.0132 17.0732 0.7743 0.0015 0.0330 0.0011
4 0.0001 265.4613 0.2193 0.9984 0.0175 0.998910장 - 공선형 데이터의 처리
data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")
head(data_9.5) YEAR IMPORT DOPROD STOCK CONSUM
1 49 15.9 149.3 4.2 108.1
2 50 16.4 161.2 4.1 114.8
3 51 19.0 171.5 3.1 123.2
4 52 19.1 175.5 3.1 126.9
5 53 18.8 180.8 1.1 132.1
6 54 20.4 190.7 2.2 137.7# 회귀분석(2) : 데이터 1949~1959
data_use<-subset(data_9.5,YEAR<=59)
res<-lm(IMPORT~.,data_use[-1]) #1열 제거
summary(res)
Call:
lm(formula = IMPORT ~ ., data = data_use[-1])
Residuals:
Min 1Q Median 3Q Max
-0.52367 -0.38953 0.05424 0.22644 0.78313
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.12799 1.21216 -8.355 6.9e-05 ***
DOPROD -0.05140 0.07028 -0.731 0.488344
STOCK 0.58695 0.09462 6.203 0.000444 ***
CONSUM 0.28685 0.10221 2.807 0.026277 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4889 on 7 degrees of freedom
Multiple R-squared: 0.9919, Adjusted R-squared: 0.9884
F-statistic: 285.6 on 3 and 7 DF, p-value: 1.112e-07head(data_use[-c(1:2)]) DOPROD STOCK CONSUM
1 149.3 4.2 108.1
2 161.2 4.1 114.8
3 171.5 3.1 123.2
4 175.5 3.1 126.9
5 180.8 1.1 132.1
6 190.7 2.2 137.7### 주성분(principle component)
pc<-prcomp(data_use[-c(1:2)],scale.=T)
pc$rotation PC1 PC2 PC3
DOPROD 0.70633041 -0.03568867 -0.706982083
STOCK 0.04350059 0.99902908 -0.006970795
CONSUM 0.70654444 -0.02583046 0.707197102pc$x #변화된 새로운 변수가 저장된곳 PC1 PC2 PC3
1 -2.1258872 0.63865815 -0.02072230
2 -1.6189273 0.55553922 -0.07111317
3 -1.1151675 -0.07297970 -0.02173008
4 -0.8942966 -0.08236998 0.01081318
5 -0.6442081 -1.30668523 0.07258248
6 -0.1903514 -0.65914745 0.02655252
7 0.3596219 -0.74367447 0.04278124
8 0.9718018 1.35405877 0.06286252
9 1.5593159 0.96404558 0.02357446
10 1.7669951 1.01521706 -0.04498818
11 1.9311034 -1.66266195 -0.08061267### 주성분회귀 - 원래데이터를 주성분분석을 통해 새로운 데이터를 생성 이를 가지고 회귀
df<-data.frame(IMPORT = scale(data_use$IMPORT),
pc$x)
df IMPORT PC1 PC2 PC3
1 -1.3185185 -2.1258872 0.63865815 -0.02072230
2 -1.2084753 -1.6189273 0.55553922 -0.07111317
3 -0.6362502 -1.1151675 -0.07297970 -0.02173008
4 -0.6142416 -0.8942966 -0.08236998 0.01081318
5 -0.6802675 -0.6442081 -1.30668523 0.07258248
6 -0.3281290 -0.1903514 -0.65914745 0.02655252
7 0.1780700 0.3596219 -0.74367447 0.04278124
8 1.0143989 0.9718018 1.35405877 0.06286252
9 1.3665374 1.5593159 0.96404558 0.02357446
10 1.2564942 1.7669951 1.01521706 -0.04498818
11 0.9703816 1.9311034 -1.66266195 -0.08061267res<-lm(IMPORT~.,data=df)
summary(res) #유의한 1,2번째 계수들만 사용하고 3번째 것은 없이 해도 되겠다..
Call:
lm(formula = IMPORT ~ ., data = df)
Residuals:
Min 1Q Median 3Q Max
-0.11525 -0.08573 0.01194 0.04984 0.17236
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.900e-16 3.244e-02 0.000 1.000000
PC1 6.900e-01 2.406e-02 28.673 1.61e-08 ***
PC2 1.913e-01 3.406e-02 5.617 0.000801 ***
PC3 1.160e+00 6.559e-01 1.768 0.120376
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1076 on 7 degrees of freedom
Multiple R-squared: 0.9919, Adjusted R-squared: 0.9884
F-statistic: 285.6 on 3 and 7 DF, p-value: 1.112e-07#중간고사 이후의것이 주가나오겠지만 앞에것도 알아야함 연습문제에서 나올것 예상11장
# 11.10 -
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")
head(data_3.3) Y X1 X2 X3 X4 X5 X6
1 43 51 30 39 61 92 45
2 63 64 51 54 63 73 47
3 71 70 68 69 76 86 48
4 61 63 45 47 54 84 35
5 81 78 56 66 71 83 47
6 43 55 49 44 54 49 34# 분산확대 인자(VIF) : 10초과 > 심각한 공산성의 문제가 있음
res<-lm(Y~.,data=data_3.3)
summary(res)
Call:
lm(formula = Y ~ ., data = data_3.3)
Residuals:
Min 1Q Median 3Q Max
-10.9418 -4.3555 0.3158 5.5425 11.5990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.78708 11.58926 0.931 0.361634
X1 0.61319 0.16098 3.809 0.000903 ***
X2 -0.07305 0.13572 -0.538 0.595594
X3 0.32033 0.16852 1.901 0.069925 .
X4 0.08173 0.22148 0.369 0.715480
X5 0.03838 0.14700 0.261 0.796334
X6 -0.21706 0.17821 -1.218 0.235577
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05olsrr::ols_vif_tol(res) #10을 넘는 것이 없기에 공산성의 문제는 없다 Variables Tolerance VIF
1 X1 0.3749447 2.667060
2 X2 0.6246520 1.600891
3 X3 0.4403263 2.271043
4 X4 0.3248624 3.078226
5 X5 0.8142600 1.228109
6 X6 0.5124025 1.951591# 수정결정계수
res_summ<-summary(res)
res_summ$adj.r.squared #0.662846[1] 0.662846# Mallow's Cp - 작을수록 좋음
# 축소모형
res_subset<-lm(Y~X1+X3,data=data_3.3)
olsrr::ols_mallows_cp(res_subset,res) [1] 1.114811# AIC - 참고만
olsrr::ols_aic(res_subset,method = "SAS")[1] 118.0024# BIC
olsrr::ols_sbic(res_subset,res)[1] 121.0938### 전진적 선택법 - AIC
res_step<-step(res,direction = "forward")Start: AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6### 후진적 제거법 -AIC
res_step<-step(res,direction = "backward")Start: AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6
Df Sum of Sq RSS AIC
- X5 1 3.41 1152.4 121.45
- X4 1 6.80 1155.8 121.54
- X2 1 14.47 1163.5 121.74
- X6 1 74.11 1223.1 123.24
<none> 1149.0 123.36
- X3 1 180.50 1329.5 125.74
- X1 1 724.80 1873.8 136.04
Step: AIC=121.45
Y ~ X1 + X2 + X3 + X4 + X6
Df Sum of Sq RSS AIC
- X4 1 10.61 1163.0 119.73
- X2 1 14.16 1166.6 119.82
- X6 1 71.27 1223.7 121.25
<none> 1152.4 121.45
- X3 1 177.74 1330.1 123.75
- X1 1 724.70 1877.1 134.09
Step: AIC=119.73
Y ~ X1 + X2 + X3 + X6
Df Sum of Sq RSS AIC
- X2 1 16.10 1179.1 118.14
- X6 1 61.60 1224.6 119.28
<none> 1163.0 119.73
- X3 1 197.03 1360.0 122.42
- X1 1 1165.94 2328.9 138.56
Step: AIC=118.14
Y ~ X1 + X3 + X6
Df Sum of Sq RSS AIC
- X6 1 75.54 1254.7 118.00
<none> 1179.1 118.14
- X3 1 186.12 1365.2 120.54
- X1 1 1259.91 2439.0 137.94
Step: AIC=118
Y ~ X1 + X3
Df Sum of Sq RSS AIC
<none> 1254.7 118.00
- X3 1 114.73 1369.4 118.63
- X1 1 1370.91 2625.6 138.16### 단계적 방법
res_step<-step(res) #최종적으로는 X1+X3이 선택됨 이거 사용하는 것이 좋을듯 Start: AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6
Df Sum of Sq RSS AIC
- X5 1 3.41 1152.4 121.45
- X4 1 6.80 1155.8 121.54
- X2 1 14.47 1163.5 121.74
- X6 1 74.11 1223.1 123.24
<none> 1149.0 123.36
- X3 1 180.50 1329.5 125.74
- X1 1 724.80 1873.8 136.04
Step: AIC=121.45
Y ~ X1 + X2 + X3 + X4 + X6
Df Sum of Sq RSS AIC
- X4 1 10.61 1163.0 119.73
- X2 1 14.16 1166.6 119.82
- X6 1 71.27 1223.7 121.25
<none> 1152.4 121.45
- X3 1 177.74 1330.1 123.75
- X1 1 724.70 1877.1 134.09
Step: AIC=119.73
Y ~ X1 + X2 + X3 + X6
Df Sum of Sq RSS AIC
- X2 1 16.10 1179.1 118.14
- X6 1 61.60 1224.6 119.28
<none> 1163.0 119.73
- X3 1 197.03 1360.0 122.42
- X1 1 1165.94 2328.9 138.56
Step: AIC=118.14
Y ~ X1 + X3 + X6
Df Sum of Sq RSS AIC
- X6 1 75.54 1254.7 118.00
<none> 1179.1 118.14
- X3 1 186.12 1365.2 120.54
- X1 1 1259.91 2439.0 137.94
Step: AIC=118
Y ~ X1 + X3
Df Sum of Sq RSS AIC
<none> 1254.7 118.00
- X3 1 114.73 1369.4 118.63
- X1 1 1370.91 2625.6 138.16summary(res_step) #X3가 의미 없다고 나와도 아님 검증된것임
Call:
lm(formula = Y ~ X1 + X3, data = data_3.3)
Residuals:
Min 1Q Median 3Q Max
-11.5568 -5.7331 0.6701 6.5341 10.3610
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.8709 7.0612 1.398 0.174
X1 0.6435 0.1185 5.432 9.57e-06 ***
X3 0.2112 0.1344 1.571 0.128
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.817 on 27 degrees of freedom
Multiple R-squared: 0.708, Adjusted R-squared: 0.6864
F-statistic: 32.74 on 2 and 27 DF, p-value: 6.058e-080613 - regression analysis
연습문제 11.5
data_use<-read.table("All_Data/p329.txt",header=T,sep="\t")
head(data_use) X1 X2 X3 X4 X5 X6 X7 X8 X9 Y
1 4.918 1 3.472 0.998 1 7 4 42 0 25.9
2 5.021 1 3.531 1.500 2 7 4 62 0 29.5
3 4.543 1 2.275 1.175 1 6 3 40 0 27.9
4 4.557 1 4.050 1.232 1 6 3 54 0 25.9
5 5.060 1 4.455 1.121 1 6 3 42 0 29.9
6 3.891 1 4.455 0.988 1 6 3 56 0 29.9### 데이터 탐색 - 자료형
sapply(data_use,class) X1 X2 X3 X4 X5 X6 X7 X8
"numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
X9 Y
"numeric" "numeric" ### 데이터 탐색 - 기초 통계량
summary(data_use) #결측값의 여부를 확인 X1 X2 X3 X4
Min. :3.891 Min. :1.000 Min. :2.275 Min. :0.975
1st Qu.:5.058 1st Qu.:1.000 1st Qu.:4.855 1st Qu.:1.161
Median :5.974 Median :1.000 Median :5.685 Median :1.432
Mean :6.405 Mean :1.167 Mean :6.033 Mean :1.384
3rd Qu.:7.873 3rd Qu.:1.500 3rd Qu.:7.158 3rd Qu.:1.577
Max. :9.142 Max. :1.500 Max. :9.890 Max. :1.831
X5 X6 X7 X8 X9
Min. :0.000 Min. :5.0 Min. :2.000 Min. : 3.00 Min. :0.00
1st Qu.:1.000 1st Qu.:6.0 1st Qu.:3.000 1st Qu.:30.00 1st Qu.:0.00
Median :1.000 Median :6.0 Median :3.000 Median :40.00 Median :0.00
Mean :1.312 Mean :6.5 Mean :3.167 Mean :37.46 Mean :0.25
3rd Qu.:2.000 3rd Qu.:7.0 3rd Qu.:3.250 3rd Qu.:48.50 3rd Qu.:0.25
Max. :2.000 Max. :8.0 Max. :4.000 Max. :62.00 Max. :1.00
Y
Min. :25.90
1st Qu.:29.90
Median :33.70
Mean :34.63
3rd Qu.:38.15
Max. :45.80 ### 데이터 탐색 - 히스토그램 & 상자그림
hist(data_use$Y,main="histogram of data_use$Y")
boxplot(data_use$Y)
### 데이터 탐색 - 산점도 행렬 & 상관계수
plot(data_use)
cor(data_use) X1 X2 X3 X4 X5 X6 X7
X1 1.0000000 0.6512669 0.6892117 0.7342737 0.45855650 0.6406157 0.3671126
X2 0.6512669 1.0000000 0.4129558 0.7285916 0.22402204 0.5103104 0.4264014
X3 0.6892117 0.4129558 1.0000000 0.5715520 0.20466375 0.3921244 0.1516093
X4 0.7342737 0.7285916 0.5715520 1.0000000 0.35888351 0.6788606 0.5743353
X5 0.4585565 0.2240220 0.2046638 0.3588835 1.00000000 0.5893871 0.5412988
X6 0.6406157 0.5103104 0.3921244 0.6788606 0.58938707 1.0000000 0.8703883
X7 0.3671126 0.4264014 0.1516093 0.5743353 0.54129880 0.8703883 1.0000000
X8 -0.4371012 -0.1007485 -0.3527514 -0.1390869 -0.02016883 0.1242663 0.3135114
X9 0.1466825 0.2041241 0.3059946 0.1065612 0.10161846 0.2222222 0.0000000
Y 0.8739117 0.7097771 0.6476364 0.7077656 0.46146792 0.5284436 0.2815200
X8 X9 Y
X1 -0.43710116 0.1466825 0.8739117
X2 -0.10074847 0.2041241 0.7097771
X3 -0.35275139 0.3059946 0.6476364
X4 -0.13908686 0.1065612 0.7077656
X5 -0.02016883 0.1016185 0.4614679
X6 0.12426629 0.2222222 0.5284436
X7 0.31351144 0.0000000 0.2815200
X8 1.00000000 0.2257796 -0.3974034
X9 0.22577960 1.0000000 0.2668783
Y -0.39740338 0.2668783 1.0000000pairs(data_use)
panel.cor<-function(x,y){
par(usr=c(0,1,0,1))
r<-round(cor(x,y),digits = 3)
text(0.5,0.5,r,cex=1.5)
}
pairs(data_use,lower.panel = panel.cor)
### (a) 모든 변수들이 모형에 포함시킬 것인가?
cor(data_use) X1 X2 X3 X4 X5 X6 X7
X1 1.0000000 0.6512669 0.6892117 0.7342737 0.45855650 0.6406157 0.3671126
X2 0.6512669 1.0000000 0.4129558 0.7285916 0.22402204 0.5103104 0.4264014
X3 0.6892117 0.4129558 1.0000000 0.5715520 0.20466375 0.3921244 0.1516093
X4 0.7342737 0.7285916 0.5715520 1.0000000 0.35888351 0.6788606 0.5743353
X5 0.4585565 0.2240220 0.2046638 0.3588835 1.00000000 0.5893871 0.5412988
X6 0.6406157 0.5103104 0.3921244 0.6788606 0.58938707 1.0000000 0.8703883
X7 0.3671126 0.4264014 0.1516093 0.5743353 0.54129880 0.8703883 1.0000000
X8 -0.4371012 -0.1007485 -0.3527514 -0.1390869 -0.02016883 0.1242663 0.3135114
X9 0.1466825 0.2041241 0.3059946 0.1065612 0.10161846 0.2222222 0.0000000
Y 0.8739117 0.7097771 0.6476364 0.7077656 0.46146792 0.5284436 0.2815200
X8 X9 Y
X1 -0.43710116 0.1466825 0.8739117
X2 -0.10074847 0.2041241 0.7097771
X3 -0.35275139 0.3059946 0.6476364
X4 -0.13908686 0.1065612 0.7077656
X5 -0.02016883 0.1016185 0.4614679
X6 0.12426629 0.2222222 0.5284436
X7 0.31351144 0.0000000 0.2815200
X8 1.00000000 0.2257796 -0.3974034
X9 0.22577960 1.0000000 0.2668783
Y -0.39740338 0.2668783 1.0000000model_1<-lm(data_use$Y~.,data_use)
summary(model_1) #전형적인 다중공선성이 보인다
Call:
lm(formula = data_use$Y ~ ., data = data_use)
Residuals:
Min 1Q Median 3Q Max
-3.7729 -1.9801 -0.0868 1.6615 4.2618
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.31044 5.96093 2.568 0.0223 *
X1 1.95413 1.03833 1.882 0.0808 .
X2 6.84552 4.33529 1.579 0.1367
X3 0.13761 0.49436 0.278 0.7848
X4 2.78143 4.39482 0.633 0.5370
X5 2.05076 1.38457 1.481 0.1607
X6 -0.55590 2.39791 -0.232 0.8200
X7 -1.24516 3.42293 -0.364 0.7215
X8 -0.03800 0.06726 -0.565 0.5810
X9 1.70446 1.95317 0.873 0.3976
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.973 on 14 degrees of freedom
Multiple R-squared: 0.8512, Adjusted R-squared: 0.7555
F-statistic: 8.898 on 9 and 14 DF, p-value: 0.0002015#coef가 X6 0.820이므로 확인 필요
olsrr::ols_vif_tol(model_1) #X6,X7이 10에 가깝기에 제거해야하는 대상중 최우선 Variables Tolerance VIF
1 X1 0.14241176 7.021892
2 X2 0.35267392 2.835480
3 X3 0.40735454 2.454864
4 X4 0.26066568 3.836332
5 X5 0.54842184 1.823414
6 X6 0.08539078 11.710866
7 X7 0.10286111 9.721847
8 X8 0.43083904 2.321052
9 X9 0.51482060 1.942424#다중공선성이 보이기에 모든 변수를 모형에 포함 시킬수는 없다
#10보다 크면 심각한 공선성이 존재
### (b) 지방세(X1) 방의 수(X6), 건물의 나이(X8)가 판매가격(Y)을
# 설명하는데 적절하다는 의견에 동의 하는가?
model_2<-lm(Y~X1+X6+X8,data=data_use)
summary(model_2)
Call:
lm(formula = Y ~ X1 + X6 + X8, data = data_use)
Residuals:
Min 1Q Median 3Q Max
-3.7486 -2.4082 -0.3594 2.1378 6.5353
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.796013 4.971105 2.976 0.007462 **
X1 3.489464 0.729368 4.784 0.000113 ***
X6 -0.415515 1.182262 -0.351 0.728921
X8 0.004923 0.063597 0.077 0.939062
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.123 on 20 degrees of freedom
Multiple R-squared: 0.7655, Adjusted R-squared: 0.7303
F-statistic: 21.76 on 3 and 20 DF, p-value: 1.653e-06# VIF
olsrr::ols_vif_tol(model_2) Variables Tolerance VIF
1 X1 0.3184367 3.140342
2 X6 0.3875669 2.580200
3 X8 0.5317389 1.880622# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(model_2)
layout(1)
# Cook의 거리 - 1을 넘으면 확인을 해야한다
olsrr::ols_plot_cooksd_chart(model_2) 
#동의는 할 수 있으나(적절하지만) 최고의 모형이라는데는 동의 못함
### (c) 지방세 X1가 단독으로 판매가격 Y을 설명하는데 적절하다는 의견에 동의?
model_3<-lm(Y~X1,data=data_use)
summary(model_3) #adj-rsquared는 변수를 선택할때 사용
Call:
lm(formula = Y ~ X1, data = data_use)
Residuals:
Min 1Q Median 3Q Max
-3.8445 -2.3340 -0.3841 1.9689 6.3005
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.3553 2.5955 5.146 3.71e-05 ***
X1 3.3215 0.3939 8.433 2.44e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.988 on 22 degrees of freedom
Multiple R-squared: 0.7637, Adjusted R-squared: 0.753
F-statistic: 71.11 on 1 and 22 DF, p-value: 2.435e-08# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(model_3)
layout(1)
# Cook의 거리 - 1을 넘으면 확인을 해야한다
olsrr::ols_plot_cooksd_chart(model_3) 
### 적절한 모형 제시
data_use_2<-data_use[-6]
model_4<-lm(Y~.,data_use_2)
summary(model_4)
Call:
lm(formula = Y ~ ., data = data_use_2)
Residuals:
Min 1Q Median 3Q Max
-3.7340 -1.8513 -0.0154 1.5472 4.2113
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.40540 4.36026 3.304 0.00482 **
X1 1.81642 0.82431 2.204 0.04360 *
X2 7.13892 4.01353 1.779 0.09556 .
X3 0.14721 0.47683 0.309 0.76177
X4 2.73339 4.24921 0.643 0.52976
X5 2.06520 1.33883 1.543 0.14377
X7 -1.91236 1.79355 -1.066 0.30318
X8 -0.03832 0.06509 -0.589 0.56481
X9 1.48746 1.65931 0.896 0.38418
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.878 on 15 degrees of freedom
Multiple R-squared: 0.8506, Adjusted R-squared: 0.771
F-statistic: 10.68 on 8 and 15 DF, p-value: 5.936e-05# VIF
olsrr::ols_vif_tol(model_4) Variables Tolerance VIF
1 X1 0.2117072 4.723504
2 X2 0.3855308 2.593827
3 X3 0.4102334 2.437637
4 X4 0.2612467 3.827800
5 X5 0.5495336 1.819725
6 X7 0.3510102 2.848920
7 X8 0.4310175 2.320091
8 X9 0.6683173 1.496295# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(model_4)
layout(1)
# Cook의 거리 - 1을 넘으면 확인을 해야한다
olsrr::ols_plot_cooksd_chart(model_4) 
### 단계적 선택법 - AIC
model_5<-step(model_4)Start: AIC=57.45
Y ~ X1 + X2 + X3 + X4 + X5 + X7 + X8 + X9
Df Sum of Sq RSS AIC
- X3 1 0.789 125.00 55.605
- X8 1 2.870 127.08 56.001
- X4 1 3.426 127.63 56.106
- X9 1 6.654 130.86 56.706
- X7 1 9.414 133.62 57.207
<none> 124.21 57.453
- X5 1 19.702 143.91 58.987
- X2 1 26.198 150.40 60.046
- X1 1 40.207 164.41 62.184
Step: AIC=55.61
Y ~ X1 + X2 + X4 + X5 + X7 + X8 + X9
Df Sum of Sq RSS AIC
- X8 1 3.369 128.36 54.243
- X4 1 4.816 129.81 54.513
- X9 1 9.581 134.58 55.378
- X7 1 9.655 134.65 55.391
<none> 125.00 55.605
- X5 1 18.957 143.95 56.994
- X2 1 25.474 150.47 58.057
- X1 1 53.245 178.24 62.122
Step: AIC=54.24
Y ~ X1 + X2 + X4 + X5 + X7 + X9
Df Sum of Sq RSS AIC
- X4 1 5.011 133.38 53.163
- X9 1 6.543 134.91 53.437
<none> 128.36 54.243
- X5 1 20.033 148.40 55.724
- X7 1 22.819 151.18 56.170
- X2 1 24.299 152.66 56.404
- X1 1 95.134 223.50 65.552
Step: AIC=53.16
Y ~ X1 + X2 + X5 + X7 + X9
Df Sum of Sq RSS AIC
- X9 1 6.223 139.60 52.257
<none> 133.38 53.163
- X5 1 17.801 151.18 54.169
- X7 1 17.873 151.25 54.181
- X2 1 39.335 172.71 57.365
- X1 1 157.972 291.35 69.915
Step: AIC=52.26
Y ~ X1 + X2 + X5 + X7
Df Sum of Sq RSS AIC
<none> 139.60 52.257
- X5 1 20.836 160.43 53.596
- X7 1 21.669 161.27 53.720
- X2 1 47.409 187.01 57.274
- X1 1 156.606 296.20 68.312summary(model_5)
Call:
lm(formula = Y ~ X1 + X2 + X5 + X7, data = data_use_2)
Residuals:
Min 1Q Median 3Q Max
-4.5605 -2.0856 0.0238 1.8580 3.8981
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.6212 3.6725 3.709 0.001489 **
X1 2.4123 0.5225 4.617 0.000188 ***
X2 8.4589 3.3300 2.540 0.019970 *
X5 2.0604 1.2235 1.684 0.108541
X7 -2.2154 1.2901 -1.717 0.102176
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.711 on 19 degrees of freedom
Multiple R-squared: 0.8321, Adjusted R-squared: 0.7968
F-statistic: 23.54 on 4 and 19 DF, p-value: 3.866e-07# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(model_5)
layout(1)
# Cook의 거리 - 0.5보다 작으면 괜찮음 / 1보다 큰 값을 고려 / 전체적으로 봤을때 튀는 값이 있는 경우
olsrr::ols_plot_cooksd_chart(model_5) 
### 17번 제거
data_use_3<-model_5$model[-17,]
model_6<-lm(Y~.,data_use_3)
summary(model_6)
Call:
lm(formula = Y ~ ., data = data_use_3)
Residuals:
Min 1Q Median 3Q Max
-4.4577 -1.6655 0.1575 1.7978 4.1865
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.2263 3.4659 3.816 0.00127 **
X1 1.7014 0.6247 2.723 0.01394 *
X2 12.0705 3.6958 3.266 0.00429 **
X5 2.4602 1.1726 2.098 0.05028 .
X7 -2.2310 1.2152 -1.836 0.08294 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.553 on 18 degrees of freedom
Multiple R-squared: 0.8418, Adjusted R-squared: 0.8067
F-statistic: 23.95 on 4 and 18 DF, p-value: 5.316e-07# VIF
olsrr::ols_vif_tol(model_6) Variables Tolerance VIF
1 X1 0.3319061 3.012900
2 X2 0.3658998 2.732989
3 X5 0.5664655 1.765333
4 X7 0.6043771 1.654596# 회귀진단 그래프들
layout(matrix(1:4,nrow=2,byrow=T))
plot(model_6)
layout(1)
# Cook의 거리 - 1을 넘으면 확인을 해야한다
olsrr::ols_plot_cooksd_chart(model_6)